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Volume I: Technical Assessment Report 
1.0 Notification and Authorization 

Mr. Daniel Murri, NASA Technical Fellow for Flight Mechanics, requested that the NASA 
Engineering and Safety Center (NESC) support an inquiry from Mr. Harold Bell, Director 
Advanced Planning and Analysis Division, NASA Office of Chief Engineer, to establish the 
Simulation Framework for Rapid Entry, Descent, and Landing (EDL) Analysis assessment. The 
principal focus of the proposed activity was to develop a simulation framework and a set of 
validated and documented subsystem models and scripts to allow rapid evaluation of EDL 
characteristics in systems analysis studies, preliminary design, mission development and 
execution, and time-critical assessments. 

The Initial Evaluation for this assessment was approved by the NESC Review Board (NRB) on 
March 12, 2009. Mr. Daniel Murri was selected to lead this assessment. The Assessment Plan 
was presented and approved by the NRB on March 26, 2009. The Phase 1 Stakeholder Outbrief 
and final report were approved by the NRB on November 5 and December 17, 2009, 
respectively. 

The Phase 2 Assessment Plan was approved by the NRB on December 17, 2009. The Phase 2 
final report was presented to the NRB for approval on February 3, 2011. 

The key stakeholders for this activity are Mr. Harold Bell; Mr. Thomas Zang, EDL-Systems 
Analysis (EDL-SA) Team Lead; and the NESC. 
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data extracted from Program/Project documents, contractor reports, and open literature, and/or 
generated from independently conducted tests, analysis, and inspections. 
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4.0 Executive Summary 

The NASA Engineering and Safety Center (NESC) was requested to establish the Simulation 
Framework for Rapid Entry, Descent, and Landing (EDL) Analysis assessment, which involved 
development of an enhanced simulation architecture using the Program to Optimize Simulated 
Trajectories II (POST2) 1 simulation tool. The assessment was requested to enhance the 
capability of the Agency to provide rapid evaluation of EDL characteristics in systems analysis 
studies, preliminary design, mission development and execution, and time-critical assessments. 
Many of the new simulation framework capabilities were developed to support the Agency EDL- 
Systems Analysis (SA) team that is conducting studies of the technologies and architectures that 
are required to enable human and higher mass robotic missions to Mars. 

The EDL-SA studies have been conducted using POST2, but with a variety of ad-hoc and 
undocumented sub-system models (e.g., mass, aerodynamics, atmosphere, guidance, control) and 
scripts. During Phase 1 of the assessment [ref.l], the NESC team developed a simulation 
framework and a set of validated and documented sub-system models (including mass models, a 
pseudo bank angle controller, aerodynamic trim, aerodynamic data models, atmospheric models, 
guidance algorithms) and scripts. Phase 2 augments the previous effort with several models that 
were not included in Phase 1, which add to the simulation support of the EDL-SA team, and 
Agency and NESC rapid assessments. These models include a navigation filter for EDL, 
aerodynamic uncertainty, second order actuator model, aerocapture guidance algorithm, 
aerobraking mission design, Earth-Global Reference Atmosphere Model (GRAM) 2010 [ref. 2], 
Earth Gravity Model (EGM) 96, and a three-axis pseudo-controller. A set of input cases was 
generated for the models and is maintained with the simulation. This test suite will be used to 
confirm that future software models added to the simulation do not adversely affect the operation 
of the existing models added from this assessment. In the POST2 simulation environment, the 
simulation framework developed during this assessment is referred to as Rapid EDL Analysis 
Simulation (REDLAS). 

The observations and NESC recommendations are discussed in Section 6.0. Overall, the main 
objective to increase the Agency’s ability to rapidly evaluate EDL characteristics in systems 
analysis studies, preliminary design, mission development and execution, and time-critical 
assessments was accomplished. It is recommended that the custodians of the POST2 simulation 
framework, the Atmospheric Flight and Entry Systems Branch at Langley Research Center 
(LaRC), should retain and maintain the REDLAS models, test cases, and scripts for use by 
current and future programs. 


POST2 is an ITAR restricted code and may only be distributed to approved persons, https://post2.larc.nasa.gov/ 
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5.0 Models and Script Descriptions 

The NESC was requested to establish the Simulation Framework for Rapid EDL Analysis 
assessment, which involved development of an enhanced simulation architecture using the 
POST2 simulation tool. The simulation developed in Phase 1 included mass models, a pseudo 
controller, aerodynamic models, atmospheric models, guidance algorithms, and various support 
scripts. Phase 2 augments the previous effort with several models that were not included in 
Phase 1, which add to the simulation support of the EDL-SA team and Agency and NESC rapid 
assessments. A set of input cases was generated for the models and is included with the 
simulation. This test suite can be used to confirm that future software models added to the 
simulation do not adversely affect the operation of the existing models from this assessment. 

The models described below were developed and/or implemented by the NESC team. This 
report is outlined and discussed as follows: 


Section 5.1: 
Section 5.2: 
Section 5.3: 
Section 5.4: 
Section 5.5: 
Section 5.6: 


multimode Extended Kalman Filter (mEKF) 

Aerodynamic Uncertainty Models 

Terminal Point Control (TPC) Guidance Algorithm 

Aerobraking Mission Design Module 

Earth Environment Models (atmosphere and gravity) 

Attitude Control Models (3-axis pseudo-controller (3APC), second order actuator 
model) 


Table 5.0.1 lists the individuals responsible for the various models and sections of this report. 
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Table 5 . 0 . 1 . Responsible Individuals for Models Implemented and Report Sections 


Model Description 

Section(s) 

Responsible Individual 

mEKF 

5. 1 and Appendix A 

Robert Bishop 

POST2 sensor models 
and input/output 

5.1.1 -5.1.2 

Scott Striepe 

Aerodynamics 
Uncertainty Models 

5.2 

Mark Schoenenberger 

POST2 Aerodynamic 
Uncertainty Inputs & 
Outputs 

5.2.1 

Scott Striepe 

TPC Guidance 
Algorithms 

5.3 

Eric Queen 

Aerobraking Mission 
Design Module 

5.4 

Alicia Cianciolo 

Earth-GRAM 2010 

5.5.1 

John Aguirre/Dick Powell 

EGM96 POST2 inputs 

5.5.2 

Scott Striepe 

Three-Axis Pseudo- 
Controller 

5.6.1 

Eric Queen/Scott Striepe 

Second Order Actuator 
Model 

5.6.2 

Eric Queen 


5.1 Multi-mode EKF 

A significant area of investigation supporting the NESC’s Simulation Framework for Rapid EDL 
Analysis is the development of a mEKF for navigation. While using perfect navigation (i.e., 
perfect knowledge of the true state) with guidance systems is a reasonable first design step, 
errors from imperfect state knowledge help determine the robustness of the guidance and vehicle 
design. The mEKF navigation filter architecture for application to EDL is consistent with the 
methodology being used for flight systems that are currently under development. 

A key element of the mEKF is sensor error modeling. The objective is to construct a generalized 
information structure to represent sensors that encompass the requisite data required to compute 
and process mEKF measurements. The mEKF navigation filter architecture can be used in 
applications of EDL using a standard set of measurements including altimeters, velocimeters, 
star trackers, and inertial measurement units. The navigation architecture main computational 
unit is the mEKF algorithm employing a dual-inertial state implementation with integrated 
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position, velocity, and attitude. The filter described here is readily re-configurable to 
accommodate differing mission scenarios (e.g., lunar descent or Mars entry). 

Also, several sensor models were included in POST2 to provide mEKF measurements. These 
models create measurements from values calculated during the simulation and applying errors 
such as bias, scale factors, and random noise. The inertial measurement unit (IMU) model 
generates A -V and A -9 measurements sensed acceleration and angular velocity. The vehicle 
altitude and velocity relative to the surface are used for the altimeter and velocimeter, 
respectively. The star tracker model uses the inertial to body transformation to define attitude 
measurements. 

Details of the theory and equations used for the mEKF and its sensor error models are provided 
in Appendix A. Sensor models in POST2 to provide measurement data are given in Section 
5.1.1. Inputs and outputs associated with these models as implemented in POST2 are identified 
in Section 5.1.2. 

5.1.1 POST2 Sensor Models for mEKF 

5. 1.1.1 IMU and Gyroscope Model 

The POST2-developed IMU model is a statistically based accelerometer and gyroscope model. 
The accelerometer model takes the true body acceleration, with a random bias, noise, and scale 
factor errors to generate a measurement of accelerations that are then converted to A-V before 
being passed to mEKF for processing. The acceleration measurement model is given by: 


®acc ~ (^3x3 + SF ccc) a env + ^ acc + Vacc’ 


and 



SF 

or acc(l t 

0 

0 


"arrt.O 


*lacc(x ) 

SF„ = 

0 

SF acciy 

0 

• Kc = 

Facciy) 

• 'hr. = 

Vacc(v) 


0 

0 

SF acci:) 


^ai <(;) 


»W) 


where a acc is the accelerometer measured acceleration vector in m/s 2 ; a env is the true body sensed 
acceleration vector in m/s“; I^ is the identity matrix; SF acc is the accelerometer scale factor 
diagonal matrix; b acc is the accelerometer bias; and f/ uii is the random noise. The accelerometer 
measurement updates begin at the start of the trajectory, and continue until touchdown at a user 
input rate. The integration of the sensed acceleration across a measurement interval defines the 
A-V: 
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tk 

A - V)* = ^ dace ■ dt 
tk - 1 


(5.1-3) 


and the accumulated A-V is the sum of the A-V generated between updates to the navigation 
filter. 

Similarly, the gyroscope model takes the true body angular rates, with a random noise, bias and 
scale factor errors to generate a measurement of angular rates that are converted to A-9 before 
being passed to navigation. Similar to the accelerometer, Equations (5.1-1), (5.1-2), and (5.1-3), 
the gyroscope angular rate measurement model is: 


= ( 7 3x3 + SF g yro K,v + bgpo + Vzyro ( ^ 

where ^ is the gyroscope measured angular rate vector in rad/s; <y m is the true (environment) 
body angular rate vector in rad/s; SF is the gyroscope scale factor diagonal matrix; b gyro is the 
gyroscope bias; and rj gwo is the random noise. The integrated angular velocity, & m , is integrated 
across a measurement interval to provide the A-9 measurement: 


tk 

A — 0^)k — £ COgyro • dt 
tk - 1 


n 


(5.1-5) 


which is then accumulated (summed) until the next call to the mEKF. The gyroscope 
measurement updates also begin at the start of the trajectory and continue until touchdown at a 
user-defined rate. 


5. 1.1.2 Star Tracker Model 

The POST2-developed star tracker measurement model is a statistically based model. The star 
tracker model takes the true attitude quaternion, bias and noise to generate a star tracker 
measurement that is passed to the mEKF for processing as an external measurement. The noise 

( Tv< ) and bias ( b sl ) are used to calculate an error quaternion ( q e(b } ) : 


Qeibji) 



□ 


(5.1-6) 


where 
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P*-n*+h' 



The star tracker measured attitude quaternion ( q st ) is determined as the product of the true 
quaternion (q env ) and error quaternion ( q e(b l]} ) yielding 


5. 1.1. 3 Altimeter and Velocimeter Models 

The measured altitude (h sen sor) determined by the altimeter sensor is modeled in POST2 as: 

hsensor = htruth ( 1 + SF h ) + ( Pnoise htruth + hnoise ) + hbias (5.1-9) 

where h tru th is the true altitude, and p no ise, h n0 ise, hbias, and SFh are the noise percentage, noise 
addition, bias, and scale factor specifications, respectively. The sensor altitude measurement 
updates begin at a user defined truth altitude above the planet surface at a data rate also input by 
the user. The sensor velocity is separated into planet relative velocity vector components 
expressed in the body coordinate system. The sensor velocity for each component (v sen sor) 
independently determined by the velocimeter is modeled as: 

^sensor = Vt ru th ( 1 + SF V ) + V no i se + Vbias (5.1-10) 


where v tru th is the true relative velocity component being modeled, and v n0 i S e, Vbi as , and SF V are the 
noise, bias, and scale factor specifications, respectively. 

5.1.2 POST2 Inputs/Outputs for mEKF 

The inputs and outputs for use with the mEKF navigation filter have been integrated into the 
POST2 input/output dictionary. The inputs are identified in Table 5.1-1 and the outputs are 
identified in Table 5.1-2. 


Table 5.1-1. POST2 Inputs for the mEKF Algorithm 


Input 

Symbol 

Units / 
Type 

Stored 

Value 

Definition 

EKF_NX 

integer 

0 

Number of states in the filter 

EKF_NP 

integer 

0 

Number of states in the covariance 
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Input 

Units / 

Stored 


Symbol 

Type 

Value 

Definition 




Number of measurements 

EKF_NM 

integer 

0 


EKF_NAV_RATE 

Hz 

0.0 

Update rate of the mEKF 

EKF_NAV_DT 

s 

0.0 

Inverse of the navigation rate 

ekf_q_f_s_o. 

nd 

1.0, 

Quaternion of orientation of planet- 

EKF_Q_E_S_ 1 , 


0.0, 

fixed reference frame with respect to the 

EKF_Q_F_S_2, 


0.0, 

planet surface reference frame 

EKF_Q_F_S_3 


0.0 


EKF_R_REF_F_X, 

m 

0.0, 

Location of the surface reference frame 

EKF_R_REF_F_Y, 


0.0, 

origin in the planet fixed-frame 

EKF_R_REF_F_Z 


0.0 


EKF_OMEGA_FI_F_X, 

m 

0.0, 

Spin axis of planet or moon at epoch 

EKF_OMEGA_FI_F_Y, 


0.0, 


EKF_OMEGA_FI_F_Z 


0.0 


EKF_Q_IE_F_0, 

nd 

1.0, 

Quaternion defining orientation of the 

EKF_Q_IE_F_ 1 , 


0.0, 

inertial reference frame with respect to 

EKF_Q_IE_F_2 , 


0.0, 

the planet-fixed reference frame at 

EKF_Q_IE_F_3 


0.0 

epoch 

EKF_QIMU_iJ, 

nd 


IMU process noise matrix 

i=l,6 


0.0 


j=l,6 




EKF_QTUNE_i_i , 

nd 

0.0 

Tuning process noise matrix 

i=l ,22 




EKF_TE 

s 

0.0 

Epoch associated with Q_IE_F 

EKF_MOON_AE 

m 

1737400 

Radius of the moon for gravity 
calculations 

EKF_MOON_APO 

m 

1735970 

Polar radius of the moon 

EKF_MOON_AEQ 

m 

1738140 

Equatorial radius of the moon 

EKF_MOON_J2 

nd 

203.32E-06 

Gravity J2 term of the moon 

EKF_MOON_MU 

mV 

4902. 8E09 

Gravitational constant of the moon 

EKF_MOON_OMEGA 

rad/s 

2.6617E-06 

Rotation rate of the moon 

EKF_MARS_AE 

m 

3397200 

Radius of Mars for gravity calculations 

EKF_MARS_APO 

m 

3376200 

Polar radius of Mars 

EKF_MARS_AEQ 

m 

3396190 

Equatorial radius of Mars 

EKF_MARS_J2 

nd 

0.0020 

Gravity J2 term of Mars 

EKF_MARS_MU 

ITT 
m / s 

4.2828376383e+13 

Gravitational constant of Mars 

EKF_MARS_OMEGA 

rad/s 

7.088218e-05 

Rotation rate of Mars 

EKF_VEL_BIAS_X, 

m/s 

0.0, 

Velocimeter measurement bias in the 

EKF_VEL_BIAS_Y, 


0.0, 

body X, Y, and Z axes 

EKF_VEL_B1AS_Z 


0.0 
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Input 

Units / 

Stored 


Symbol 

Type 

Value 

Definition 

EKF_VEL_N OISE_C O V_X_ 

nd 

0.0, 

Velocimeter measurement noise 

X, 


0.0, 

covariance matrix 

EKF VEL NOISE COV X 


0.0, 


Y, 


0.0, 


EKF VEL NOISE COV X 


0.0, 


z, 


0.0, 


EKF VEL NOISE COV Y 


0.0, 


X, 


0.0, 


EKF VEL NOISE COV Y 


0.0 


Y, 




EKF VEL NOISE COV Y 




z, 




EKF VEL NOISE COV Z 




X, 




EKF VEL NOISE COV Z 




Y, 




EKF VEL NOISE COV Z 




Z 




EKF_ATTC_BIAS_X, 

m/s 

0.0, 

Attitude camera measurement bias in 

EKF_ ATTC_BIAS_Y, 


0.0, 

the body X, Y, and Z axes 

EKF_ ATTC _BIAS_Z 


0.0 


EKF_ ATTC_ 

nd 

0.0, 

Attitude camera measurement noise 

NOISE_COV_X_X, 


0.0, 

covariance matrix 

EKF_ ATTC 


0.0, 


_NOISE_COV_X_Y, 


0.0, 


EKF_ ATTC 


0.0, 


_NOISE_COV_X_Z, 


0.0, 


EKF_ ATTC 


0.0, 


_NOISE_COV_Y_X, 


0.0, 


EKF_ ATTC 


0.0 


_NOISE_COV_Y_Y, 




EKF_ ATTC 




_NOISE_COV_Y_Z, 




EKF_ ATTC 




_NOISE_COV_Z_X, 




EKF_ ATTC 




_NOISE_COV_Z_Y, 




EKF_ ATTC 




_NOISE_COV_Z_Z 




EKF_ALT_B IAS 

m 

0.0 

Spherical altimeter measurement bias 
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Input 

Symbol 

Units / 
Type 

Stored 

Value 

Definition 

EKF_ALT_N OISE_CO V 

nd 

0.0 

Spherical altimeter measurement noise 
covariance 

EKF_ACCEL_SPEC_BIAS 

m/s 

0.0 

IMU A-V bias standard deviation (1-g) 

EKF_ACCEL_SPEC_NOISE 

nd 

0.0 

IMU A-V noise covariance 

EKF_GYRO_SPEC_BIAS 

rad 

0.0 

IMU A-6 bias standard deviation (1-g) 

EKF_GYRO_SPEC_NOISE 

nd 

0.0 

IMU A-6 noise covariance 

ekf_xi_delta_err, 

m 

0.0, 

Initial navigation position error input as 

ekf_yi_delta_err, 

EKF_ZI_DELTA_ERR 


0.0, 

0.0 

distance from true inertial state 

EKF_VXI_DELT A_ERR, 

m/s 

0.0, 

Initial navigation velocity error input 

EKF_V YI_DELT A_ERR, 

ekf_vzi_delta_err 


0.0, 

0.0 

relative to true inertial state 

EKF_ERR_Q_0, 

nd 

1.0, 

Initial quaternion defining the inertial 

EKF_ERR_Q_1, 

EKF_ERR_Q_2, 

EKF_ERR_Q_0 


0.0, 

0.0, 

0.0 

navigation attitude error 
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Table 5.1-2. POST2 Outputs for the mEKF Algorithm 


Output 

Symbol 

Type/ 

Units 

Definition 

ekf_x_km_x, 

ekf_x_km_y, 

ekf_x_km_z 

m 

Propagated navigation inertial position estimate 

ekf_x_km_vx, 

ekf_x_km_vy, 

ekf_x_km_vz 

m/s 

Propagated navigation inertial velocity estimate 

ekf_x_km_q0, 

ekf_x_km_ql, 

ekf_x_km_q2, 

ekf_x_km_q3 

nd 

Propagated navigation inertial to body quaternion estimate 

ekf_P_km 

27~2 
m /s 

Propagated navigation covariance matrix 

ekf_x_kp_x, 

ekf_x_kp_y, 

ekf_x_kp_z 

m 

Navigation inertial position estimate after update 

ekf_x_kp_vx, 

ekf_x_kp_vy, 

ekf_x_kp_vz 

m/s 

Navigation inertial velocity estimate after update 

ekf_x_kp_q0, 

ekf_x_kp_ql, 

ekf_x_kp_q2, 

ekf_x_kp_q3 

nd 

Navigation inertial to body quaternion estimate after update 

ekf_P_kp 

27~2 
m /s 

Navigation covariance matrix after update 


5.2 Aerodynamic Uncertainty Model 

Part of the Phase 1 effort was to include several aerodynamic models (i.e., the 70-, 60-, and 45- 
degree sphere-cone aerodynamic database routines). However, the uncertainty values (used in 
Monte Carlo assessments) were required to be user input. This approach requires the user to 
have detailed knowledge of the uncertainty values and their application for the given database 
being used. To avoid user input errors associated with aerodynamic uncertainties, the method 
described below was developed wherein the uncertainty value and application are held within the 
aerodynamic database subroutine. The user inputs the desired standard deviation value (i.e., 
sigma value) of the uncertainty desired. The actual dispersion value and its application to the 
nominal aerodynamic coefficient values throughout the trajectory simulation are computed 
within the aerodynamic subroutine. This method is now used for the dispersion calculation 
associated with the Phase 1 aerodynamic models listed above. 

The dispersions are input in specific flight regimes. The static coefficient dispersions are defined 
in three main regions of the flight regime: free molecular, hypersonic continuum, and supersonic 
continuum. The free molecular flight regime is defined by Knudsen number values > 1000 
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(indicated by UNCI in the variable name). For Knudsen number < 1000, the hypersonic 
continuum is defined by Mach numbers between 5 and 10 ( UNC2 in variable name) and the 
supersonic continuum region is defined by Mach < 5 (_UNC3 in variable name). The 
transitional aerodynamic uncertainties are determined using linear interpolation. Interpolation 
using log base 10 of the Knudsen number is applied between free molecular and supersonic 
continuum dispersion values. Whereas, interpolation using Mach number is utilized between 
hypersonic and supersonic continuum dispersion values. The only difference for dynamic 
derivative coefficient dispersions is that the hypersonic continuum is defined by Mach numbers 
below 6, and supersonic continuum is Mach number 3 and less. 

The dispersions are determined in the DISPERSION subroutine (in the file dispersion. f in the 
redl Aero directory). This routine receives nominal aerodynamic coefficients (at respective 
database moment reference points) from the aerodynamic database routines and transforms 
moments to the center of gravity (eg). The aerodynamic coefficients are then dispersed about the 
eg using 3-sigma values [refs. 3, 4, 5] scaled by the user input (as described in Table 5.2-1). The 
dispersed aerodynamic coefficients are then transformed to the moment reference point for 
correct dispersion application by POST2 computations, which include the contribution of the 
aerodynamic forces to the aerodynamic moments. User input values of zero apply no dispersion 
and values may range from -1 to 1. Three-sigma values are multiplied by these input values (i.e., 
an input of 1 applies a +3 sigma dispersion as indicted in Table 5.2-1). 

The general variables in Table 5.2-1 are associated with the aerodynamic inputs, whereas Table 
5.2-2 are the outputs in POST2. As shown parenthetically in Table 5.2-1, similar dispersed 
values are available for the other flight regimes with the variable name changed (specifically the 
UNCI portion) as described previously. The only exception to this information is the 
CMQCWRUNC3MULT. 
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Table 5.2-1. POST2 Inputs for the Aerodynamic Dispersion Model 


Input 

Symbol 

Units 

Stored 

Value 

Definition 

CAT_UNC1_M 

ULT 

(UNC2 &UNC3) 

nd 

0.0 

Axial force coefficient multiplicative uncertainty at eg 
in free molecular flight regime. For IDISP_AERO= 1 , 
input dispersion value as fraction of 3 -sigma value; 
that is +1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of 
the 3 -sigma value. 

(Hypersonic and Supersonic flight regimes) 

CAT_UNC1_ 

ADD 

(UNC2 &UNC3) 

nd 

0.0 

Axial force coefficient additive uncertainty at eg in 
free molecular flight regime. For IDISP_AERO= 1 , 
input dispersion value as fraction of 3 -sigma value; 
that is +1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of 
the 3 -sigma value. 

(Hypersonic and Supersonic flight regimes) 

CMQ_CWR_ 

UNC1_ADD 

(UNC2 &UNC3) 

nd 

0.0 

Pitch damping and yaw moment damping due to yaw 
rate additive dispersion values at eg for free molecular 
flight regime. Input dispersion value as fraction of 3- 
sigma value; that is +1 is + 3-sigma, -1 is -3-sigma and 
0.5 is half of the 3-sigma value. Same value is applied 
to both dynamic derivatives. 

(Hypersonic and Supersonic flight regimes) 

CMQ_CWR_ 

UNC3_MULT 

nd 

0.0 

Pitch damping and yaw moment damping due to yaw 
rate multiplicative dispersion values at eg in 
supersonic flight regime. Input dispersion value as 
fraction of 3-sigma value; that is +1 is + 3-sigma, -1 is 
-3-sigma and 0.5 is half of the 3-sigma value. Same 
value is applied to both dynamic derivatives. 
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Input 


Stored 


Symbol 

Units 

Value 

Definition 




Pitching moment coefficient multiplicative uncertainty 

CMT_UNC 1_M 

nd 

0.0 

at eg in free molecular flight regime. For 

ULT 



IDISP_AERO= 1 , input dispersion value as fraction of 
3-sigma value; that is +1 is + 3-sigma, -1 is -3-sigma 
and 0.5 is half of the 3-sigma value. 

(UNC2 &UNC3) 



(Hypersonic and Supersonic flight regimes) 




Pitching moment coefficient additive uncertainty at eg 

CMT_UNC1_ 

nd 

0.0 

in free molecular flight regime. For IDISP_AERO= 1 , 

ADD 



input dispersion value as fraction of 3 -sigma value; 
that is +1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of 
the 3 -sigma value. 

(UNC2 &UNC3) 



(Hypersonic and Supersonic flight regimes) 




Normal force coefficient multiplicative uncertainty at 

CNT_UNC1_M 

nd 

0.0 

eg in free molecular flight regime. For 

ULT 



IDISP_AERO= 1 , input dispersion value as fraction of 
3-sigma value; that is +1 is + 3-sigma, -1 is -3-sigma 
and 0.5 is half of the 3-sigma value. 

(UNC2 &UNC3) 



(Hypersonic and Supersonic flight regimes) 




Normal force coefficient additive uncertainty at eg in 

CNT_UNC1_ 

nd 

0.0 

free molecular flight regime. For IDISP_AERO= 1 , 

ADD 



input dispersion value as fraction of 3 -sigma value; 
that is +1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of 
the 3 -sigma value. 

(UNC2 &UNC3) 



(Hypersonic and Supersonic flight regimes) 




Yawing moment coefficient multiplicative uncertainty 

CW_UNC1_ 

nd 

0.0 

at eg in free molecular flight regime. For 

MULT 



IDISP_AERO= 1 , input dispersion value as fraction of 
3-sigma value; that is +1 is + 3-sigma, -1 is -3-sigma 
and 0.5 is half of the 3-sigma value.. 

(UNC2 &UNC3) 



(Hypersonic and Supersonic flight regimes) 
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Input 

Symbol 

Units 

Stored 

Value 

Definition 

CW_UNC1_ 

ADD 

(UNC2 &UNC3) 

nd 

0.0 

Yawing moment coefficient additive uncertainty at eg 
in free molecular flight regime. For IDISP_AERO= 1 , 
input dispersion value as fraction of 3 -sigma value; 
that is +1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of 
the 3 -sigma value. 

(Hypersonic and Supersonic flight regimes) 

CY_UNC1_ 

MULT 

(UNC2 &UNC3) 

nd 

0.0 

Side force coefficient multiplicative uncertainty at eg 
in free molecular flight regime. For IDISP_AERO= 1 , 
input dispersion value as fraction of 3 -sigma value; 
that is +1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of 
the 3 -sigma value. 

(Hypersonic and Supersonic flight regimes) 

CY_UNC1_ 

ADD 

(UNC2 &UNC3) 

nd 

0.0 

Side force coefficient additive uncertainty at eg in free 
molecular flight regime. For IDISP_AERO=l, input 
dispersion value as fraction of 3-sigma value; that is 
+ 1 is + 3-sigma, -1 is -3-sigma and 0.5 is half of the 3- 
sigma value. 

(Hypersonic and Supersonic flight regimes) 

IDISP_AERO 

Integer 

0 

Flag to generate dispersed aerodynamics. 

=0, use nominal (un-dispersed) aerodynamic 
coefficients 

= 1 , disperse aerodynamics at eg using inputs as 
fractions of + 3-sigma dispersions. 

REDL_AEROU 
NCERTi, 
j=l ,50 

nd 

0.0 

Array of aerodynamic uncertainties applied to nominal 
coefficients. Values input for standard or total 
aerodynamic uncertainties defined in this table also 
populate this array in the appropriate location. 
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Table 5.2-2. POST2 Outputs from the Aerodynamic Dispersion Model 


Output 

Symbol 

Type/ 

Units 

Definition 

CA_NOM 

nd 

Nominal Axial force coefficient 

CAT_NOM 

nd 

Nominal total Axial force coefficient 

CLL_NOM 

nd 

Nominal rolling moment coefficient 

CLLP_NOM 

nd 

Nominal roll moment damping coefficient due to roll rate 

CLLR_NOM 

nd 

Nominal roll moment damping coefficient due to yaw rate 

CM_NOM 

nd 

Nominal pitching moment coefficient 

CMQ_NOM 

nd 

Nominal pitch moment damping coefficient due to pitch rate 

CMT_NOM 

nd 

Nominal total pitching moment coefficient 

CN_NOM 

nd 

Nominal Normal force coefficient 

CNT_NOM 

nd 

Nominal total Normal force coefficient 

CW_NOM 

nd 

Nominal yawing moment coefficient 

CWP_NOM 

nd 

Nominal yaw moment damping coefficient due to roll rate 

CWR_NOM 

nd 

Nominal yaw moment damping coefficient due to yaw rate 

CY_NOM 

nd 

Nominal Side force coefficient 

CA_DIS 

nd 

Dispersed Axial force coefficient 

CAT_DIS 

nd 

Dispersed total Axial force coefficient 

CLL_DIS 

nd 

Dispersed rolling moment coefficient 

CLLP_DIS 

nd 

Dispersed roll moment damping coefficient due to roll rate 

CLLR_DIS 

nd 

Dispersed roll moment damping coefficient due to yaw rate 

CM_DIS 

nd 

Dispersed pitching moment coefficient 

CMQ_DIS 

nd 

Dispersed pitch moment damping coefficient due to pitch rate 

CMT_DIS 

nd 

Dispersed total pitching moment coefficient 

CN_DIS 

nd 

Dispersed Normal force coefficient 

CNT_DIS 

nd 

Dispersed total Normal force coefficient 

CW_DIS 

nd 

Dispersed yawing moment coefficient 

CWP_DIS 

nd 

Dispersed yaw moment damping coefficient due to roll rate 

CWR_DIS 

nd 

Dispersed yaw moment damping coefficient due to yaw rate 

CY_DIS 

nd 

Dispersed Side force coefficient 

KNUDSEN 

nd 

Knudsen Number 


5.3 Terminal Point Control Guidance Algorithm 

A key element in developing rapid vehicle simulations is the vehicle guidance, navigation, and 
control (GN&C). The vehicle guidance (i.e., the "driver") takes input from the navigation 
system and user input targeting information to send signals to the flight control system that will 
guide the vehicle to its destination while remaining within the operating constraints. Several 
guidance algorithms were included in Phase 1. For this phase, the TPC guidance algorithm for 
aerocapture was added. 
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This guidance algorithm was implemented in order to examine vehicles that enter the atmosphere 
and aerocapture into orbit prior to EDL. The TPC guidance formulation was derived using the 
calculus of variations, which was the same methodology used to develop the Apollo Earth entry 
terminal phase guidance. A reference trajectory is used to determine the sensitivities of the final 
vehicle state to changes in the control at any point along the trajectory. For the TPC guidance, 
the control variable is the vehicle bank angle. The derivation differs from that of the Apollo 
terminal phase guidance in two significant ways. First, the boundary conditions for aerocapture 
are different from those for landing. Second, the aerocapture trajectory is not monotonic in 
altitude, so the simplification allowed by parameterization on altitude is not possible for 
aerocapture. The TPC guidance provides bank angle commands to reach a target orbit after a 
single atmospheric pass. 

As in the entry guidance, the bank angle magnitude controls the in-plane lift. For aerocapture, 
the in-plane lift targets to the desired apoapsis while maintaining periapsis as high as possible. 
The bank angle sign is used for out-of-plane control of either inclination or wedge angle. The 
TPC algorithm uses a lateral corridor that is indexed on energy. Whenever the wedge or 
inclination angle error is outside the corridor, the bank angle sign is negated to drive it back into 
the corridor. There is no general procedure to determine the shape of the corridor. It is tuned to 
meet the needs of each individual mission. The TPC guidance algorithm theoretical 
development with equations is provided in references 6 and 7. 

The process of tuning the TPC guidance starts with performing a reference 3-degree of freedom 
(DoF) trajectory. The reference trajectory is run using the input file PICKBANKNOM.INP. 

This input file can be set either to pick the required constant bank angle for a given entry flight 
path angle, or to pick the required entry flight path angle for a given bank angle profile. The 
latter allows a variable bank angle profile, typically a bank angle is used that is more lift up 
initially, more lift down at the end, and varies linearly with energy. The input file should be a 
nominal aerocapture trajectory except that the lift is in-plane and unguided (i.e., only an open 
loop bank angle profile is used). The input files PICKBANKSTP.INP and 
PICKBANKSHW.INP are the same except that one is steeper at entry, whereas the other is 
shallower. The amount of steepness or shallowness depends on the bank profile used. Targeting 
near, but not at the edge of the flyable corridor works best. Note that aerocapture trajectories can 
be sensitive. An initial bank angle estimate that results in the hyperbolic approach trajectory 
capturing into an elliptical orbit is required. A bank angle yielding in a flight path angle that is 
too steep at atmospheric interface results in the vehicle impacting into the planet’s surface, and 
thus an acceptable derivative can not be generated. If the flight path angle is too shallow, then 
the vehicle skips out (i.e., does not achieve elliptical orbit) and no acceptable derivative can be 
generated. Fortunately, first order calculations can be used to converge on initial parameters. 

After the reference trajectories are performed, form the .mat files with the extensions “_nom”, 

“ stp” and “_shw” for the nominal, steep and shallow trajectory cases, respectively. Then run 
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the script file “multimaster.m” in MATLAB™. For each of the three references, multimaster 
performs several functions: 

• Retrieves the relevant information from the .mat file and saves it to a .dat file 
(pbnom.dat, pbstp.dat and pbshw.dat). 

The following variables are needed in the print block: TIME, ALTITO, 
GAMIRV, VELR, RDOT, DRAGW, DENS, BNKANG, AEROCAP_DV, 
ENERGY, GCRAD, GXI, GYI, GZI, DECLN, LONG, ENRGDT, MALTA, 
WEIGHT, VELI, and CDW; 

• Identifies the terminal costate boundary conditions based on the exit states and the target 
orbit; 

• Integrates the costates backwards from the final condition and calculates the gains along 
the trajectory; 

• Reduces the number of data points to minimize storage requirements; 

• Plots the resulting data to check that the reduced data represents the original; and 

• Writes the gains in the form of POST2 tables, including comments. 

The gain tables are saved in the files TPC_gains_nom.dat, TPC_gains_shw.dat, and 
TPC_gains_stp.dat. The multimaster script calls: nomgenmulti.m, mgacdmulti.m, 
data_reduct.m, verify_plots.m, writepostinput_nom.m, writepostinput _stp.m, and writepostinput 
_shw.m. The script mgacdmulti.m is where the boundary conditions are calculated, the costates 
are integrated, and the gains are calculated. To run POST2 with the guidance in the loop, the 
input deck needs to include the three gain table files and the guidance setup file, tpc_setup.dat, 
and to have the guidance mode set to 15 (IGUID(14) = 15). Note that the parameters most 
adjusted to achieve a successful aerocapture are TPC_OVERKD, TPC_OVERKRD, and 
TPC_NSAMP. TPC_OVERKV is often set to zero. 

POST2 inputs and outputs for the TPC guidance algorithm are provided in Tables 5.3-1 and 
5.3-2, respectively. Note that all of the variables in the TPC guidance structure are prefixed with 
“tpc_”. Many of the input parameters are generated via MATLAB™ scripts and the nominal 
trajectory. Additionally, all of the gain tables are generated using MATLAB™ scripts. These 
gain table and input value files generated via the process must be included in the POST2 input 
file for the TPC guidance to function properly. 
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Table 5.3-1. POST2 Inputs for the TPC Guidance Algorithm 


Input 

Symbol 

Units 

Stored 

Value 

Definition 

IGUID(14) 

integer 

0 

The TPC selection flag. 
= 15, Use TPC guidance 

TPC_REQ 

m 

0.0 

Planet equatorial radius used in TPC 

TPC_RPL 

m 

0.0 

Planet polar radius used in TPC 

TPC_WIE 

rad/s 

0.0 

Planet rotation rate used in TPC 

TPC_XMUE 

3/2 

m / s 

0.0 

Planet gravitational constant used in TPC 

TPC_DT 

s 

0.0 

TPC guidance call update time 

TPC_GSTART, 

TPC_GSTOP 

m/s 2 

0.0, 

0.0 

Acceleration triggers to begin and end atmospheric 
phase of TPC guidance 

TPC_TAU1, TPC_TAU2 

s 

0.0 

TPC drag filter time constants 

TPC_CDNOM 

nd 

0.0 

Nominal vehicle drag coefficient used by TPC 

TPC_SNOM 

1 

m 

0.0 

Nominal vehicle reference area used by TPC 

TPC_DENSFACT 

nd 

0.0 

TPC over-density factor applied to the atmospheric 
density 

TPC_ 

SCALEHEIGHT_NOM 

km 

0.0 

Atmospheric scale height used to initialize the 
atmospheric density estimator in TPC 

TPC_OVERKD 

nd 

0.0 

TPC overcontrol gain on drag acceleration 

TPC_OVERKRD 

nd 

0.0 

TPC overcontrol gain on radius rate of change 

TPC_OVERKV 

nd 

0.0 

TPC overcontrol gain on velocity 

TPC_DELT A_APO 

m 

0.0 

Difference from the target apopasis radius 
(TPC_RA_TARG) at which aeroshell is jettisoned 
(e.g., if =5000, aeroshell will jettison 5 km above 
target apoapsis) 

TPC_R A_T ARG 

m 

0.0 

The target apopasis radius at which aeroshell is 
jettisoned in TPC 

TPC_RHOO 



Atmospheric surface density used to initialize the 
atmospheric density estimator in TPC 

TPC_NSAMPL 

integer 

0 

Number of acceleration measurements used by the 
atmospheric density estimator in TPC 

TPCJAXIALFLAG 

integer 

0 

Flag to determine which accelerations to use to 
estimate density in TPC 
= 2, use drag accelerations 

TPC_MIN_RHO0, 

TPC_MAX_RHOO 


0.0, 

0.0 

Fower and upper limits on the “sea-level” atmospheric 
density in the estimator in TPC 
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Input 

Symbol 

Units 

Stored 

Value 

Definition 

TPC_MIN_SCH, 

km 

0.0, 

Lower and upper limits on the atmospheric scale 

TPC_MAX_SCH 


0.0 

height in the estimator in TPC. 

TPC_IFLG 

integer 

0 

TPC mode flag. If -2, then the vehicle is forced to 
follow input reference bank profile. 

TPC_ 

BOTTOMETRIGGER 

kg nr/s 2 

0.0 

Energy trigger to change the roll reversal limits from 
TPC_B OTT OMHIE to TPC_B OTTOMLOE in TPC. 

TPC_BOTTOMLOE, 

degree 

0.0, 

The low-energy and high-energy limits on bank 

TPC_B OTT OMHIE 


0.0 

reversal direction in TPC. If energy is more than 
TPC_BOTTOMETRIGGER and 
TPC_BOTTOMHIE=80, then reversals commanded 
when bank is less than 80 will go over the top, 
reversals commanded when bank is greater than 80 
will go underneath. Set to 90 will give shortest 
distance. 

TPC_IDBG 

integer 

0 

TPC debug flag. Set =1 for additional output and 
generate an external file (fort. 68) with more output. 

TPC_LATFLAG 

integer 

0 

TPC flag to select the lateral logic variable. 
= 0, inclination 
=1, wedge angle 

tpc_banklim 

degree 

0.0 

Limits closest command to lift up or lift down (e.g., if 
tpc_banklim=10 no commands will be issued between 
10 and -10, or above 170 or below -170). 

TPC_N OMONL YFL AG 

integer 

0 

TPC flag to use only the nominal reference profile. 

TPC_TOPINCX, 

kg m 2 /s 2 

0.0, 

TPC tables used to control roll reversals. The X tables 

TPCJTOPINCY, 

0.0, 

are energy values and the Y tables are wedge angle (or 

TPC_BOTINCX, 


0.0, 

inclination) values. Whenever the vehicle exceeds the 

TPC_BOTINCY 


0.0 

interpolated Y values, a reversal is commanded. The 
tables currently have length 6. 

TPC_XMASS 


0.0 

Mass of the vehicle used in TPC 

TPC_NRDTl_NOM, 

integer 

0.0, 

Array lengths for the gains in each of the sets of tables 

TPC_NRDT1_STP, 

TPC_NRDT1_SHW 


0.0, 

0.0 

in TPC 
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Table 5.3-2. POST2 Outputs for the TPC Guidance Algorithm 


Output 

Symbol 

Type/ 

Units 

Definition 

TPC_BANKC 

degree 

Commanded Bank Angle 

TPC_JETTISO 

N_ON 

nd 

Flag indicating heatshield is to be jettisoned. 

TPC_DIR 

nd 

Commanded bank maneuver direction 
= -2, bank through 1 80 degree (underneath) 
= - 1 , bank left 
= 0, go shortest distance 
= 1 , bank right 

= 2, bank through 0 degree (over) 


5.4 Aerobraking Mission Design Module 

The program contains basic algorithms required to perform an aerobraking mission design 
simulation. It is a variation of the mission design software used for mission planning and 
aerobraking operations for Mars Odyssey and Reconnaissance Orbiter. The code allows for 
specific walk-in and operational maneuver strategies. Walk-in is the phase that transitions the 
spacecraft from its exo-atmospheric orbit to aerobraking without exceeding any of the spacecraft 
constraints (e.g., heat rate, temperature, dynamic pressure). This module is capable of 
monitoring eclipse times, lifetime constraints, and includes several options for aerobraking 
corridor control. 

The input requires the use of repeating roving events for critical events such as atmosphere entry, 
periapsis, atmosphere exit, and apoapsis. Aerobraking missions with more than 1000 orbits will 
require multiple sets of critical events in the input file. The aerobraking mission design option is 
operating when npc(9) = 3 and ab_flag=l. 

The general setup and procedure for using the aerobraking mission design module involves 
input files and scripts. Samples are provided in a redlasim_aerobraking directory. The contents 
of that directory include the following files: 

TEST.INP is a standard input used with an aerobraking case. The event structure can be 
changed, but AB_EVENT_FLAG input must adhere to the following values: 
AB_EVENT_FLAG = 30 in the atmosphere entry event; 

AB_EVENT_FLAG = 50 in the periapsis event; 

AB_EVENT_FLAG = 70 in the atmosphere exit event; and 
AB_EVENT_FLAG = 90 in the apoapsis event. 
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SPAWNP2.pl is called from PREDICT_DV.f. This script submits the POST2 spawned job 
that runs while the original aerobraking simulation is suspended. Once completed, data is 
placed into appropriate files by the script for the original run to continue. 

The RUNOUT4.PL script converts the variables saved in save_runout_variables.cpp to a 
MATLAB™ format once an input case has completed. This script creates the 
TEST_APO.mat and TEST_PER.mat files containing apoapsis and periapsis data, 
respectively. 

The RUN_S_POST.PL script submits the jobs to the cluster for execution. Currently, only 
cluster execution is permitted (no local machine runs) due to the spawn case setup is for 
cluster use only at this time. 

PLOT_AB.m is a sample plotting routine is also included for use with MATLAB™. 

This script loads the TEST.mat, TEST_PER.mat, and TEST_APO.mat to plot the results of 
the aerobraking mission. 

Test cases were developed which employ this aerobraking mission design module. The P0ST2 
inputs, tables, and outputs for the aerobraking mission design module is provided in Tables 
5.4-1, 5.4-2, and 5.4-3, respectively. 


Table 5.4-1. POST2 Inputs for the Aerobraking Mission Design Module 


Input 

Symbol 

Units 

Stored 

Value 

Definition 

npc(9) 

integer 

0 

Propulsion type selection flag 
= 0, no thrust 

= 1,2, Rocket or ramjet engine 

= 3, Instantaneous delta V option also Calculate 
Aerobraking DV 

= 4, Instantaneous delta V addition using current 
weight of propellant and specific impulse 

ABENDVALUE 

decimal 

0.0 

Value of AB END CRITR MODE to end simulation 
(i.e., ifAB END CRITR MODE = AB MALTA, 
AB END VALUE = 450) 

ABENDFLAG 

integer 

0 

Flag to end aerobraking simulation 
= 1 ends sim 

AB EVENT FLG 

integer 

0 

Aerobraking event flag 

ABFLAG 

integer 

0 

Aerobraking Mission Design module activation flag. 
=1, module active 
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ABDVMODE 

integer 

0 

AV magnitude selection flag 

0 = NO MANEUVER 

1 = TABLE LOOKUP 

2 = PREDICT 

SPAWNENDVALUE 

decimal 

0.0 

Value of SPAWN END CRITR MODE to end 
spawned case. 

ABPREDICTMODE 

integer 

0 

Type of AB corridor 

1 = ALT CORRIDOR 

2 = DENS CORRIDOR 

3 = HTRT CORRIDOR 

4 = DYNP CORRIDOR 

SPAWN END CRITR 
MODE 

integer 

0 

End of simulation criteria 

1 = ALT 

2 = ORB 

3 = DAY 

SPAWNP2PATH 

character 

*132 

Path to the SPAWNP2 to be used in the spawn cases. 
Note that the path to executable must be changed 
inside SPAWNP2 file. 

LAST WALK1N ORB 

integer 

0 

Value of the last walk in orbit. 

LAST WALKIN2 ORB 

integer 

0 

Value of last walk in orbit after solar conjunction. 


Table 5.4-2. POST2 Table Inputs for the Aerobraking Mission Design Module 


Input 

Symbol 

Units 

Stored 

Value 

Definition 

LOWERCORRT 

decimal 

0.0 

Aerobraking lower corridor limit. 

UPPERCORRT 

decimal 

Stored 

Value 

Aerobraking upper corridor limit. 

DAY SB4ABMT 

decimal 

0.0 

Days between aerobraking maneuvers. 

CORRTARGT 

decimal 

0.0 

Aerobraking corridor target (i.e., 0.5 targets the 
middle of the corridor). 

AWALKINDVT 

decimal 

0.0 

Table of walk-in maneuvers as a function of orbit 
number. 

WALKIN2DVT 

decimal 

0.0 

Table of second walk-in maneuvers as a function of 
orbit number. 
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Table 5.4-3. POST2 Outputs for the Aerobraking Mission Design Module 


Output 

Symbol 

Type/ 

Units 

Definition 

ORBITNUM 

nd 

orbit revolution counter from first walkin. Orbit numbers start 
at periapsis 

ORBIT2NUM 

nd 

Orbit revolution counter from second walkin 

ABDAYNUM 

days 

24-hr period counter (days of aerobraking) 

DOWN MANEUVER 
DV 

m/s 

down maneuver indicator, equal to DV magnitude (abs value) 

UP MANEUVER 
DV 

m/s 

up maneuver indicator, equal to DV magnitude (abs value) 

SUM MANEUVER D 
V 

m/s 

Cumulative maneuver DV 

AB MANEUVER 
COUNTER 

nd 

Orbit maneuver counter 

SIMSTARTTIME 

s 

Simulation start time 

LOWERCORR 


Value of the lower aerobraking corridor. Type of corridor and 
units depend on AB PREDICT MODE 

UPPERCORR 


Value of the upper aerobraking corridor. Type of corridor and 
units depend on AB PREDICT MODE 

CORRTARG 


Target in the corridor. Type of corridor and units depend on 
AB PREDICT MODE 

ABMAXDENS 

kg/m 3 

Maximum density on a particular aerobraking pass 

ABMAXHTRT 

w/cm 2 

Maximum heat rate indicator (1/2 rho vel A 3) on a particular 
aerobraking pass 

ABMAXDYNP 

Pa 

Maximum dynamic pressure on a particular aerobraking pass 

ABENTRYLS 

degree 

Local Solar Longitude of a particular aerobraking pass entry 


5.5 Earth Environment Models 

5.5.1 2010 Earth-Global Reference Atmosphere Model 

Earth-GRAM 2010 [ref. 2] uses an empirical database to provide atmospheric quantities (e.g., 
density, temperature, pressure, winds, and constituent concentrations) from the Earth’s surface to 
orbital altitudes based on geographic location and time of year. This version is an update to the 
previously released Earth-GRAM 2007. Specific changes for this latest release include an 
updated lower atmosphere database (Global Upper Air Climatic Atlas replaced by monthly 
global climatology by the National Centers for Environmental Prediction, or NCEP), revised 
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boundary layer model, updated thermosophere models, scaling parameters for standard 
deviations modified, auxiliary profiles option, moisture corrections, and additional output 
parameters. The Earth-GRAM 2010 models have been implemented into POST2 and are 
accessible via the standard input (npc(5)=19). The POST2 inputs, tables, and outputs for Earth- 
GRAM 2010 are listed in Tables 5.5-1, 5.5-2, and 5.5-3, respectively. 


Table 5.5-1. POST2 Inputs for the Earth-GRAM 2010 Model 


Input 

Symbol 

Type/ 

Units 

Stored 

Value 

Definition 

NPC(5) 

integer 

2 

The atmosphere model selection flag. 
=19, Earth GRAM 2010 model 

NPC(6) 

integer 

0 

Wind calculation flag 

= 4 Use GRAM determination of North-South, East- 
West, and Vertical winds 

AP 

decimal 

0.0 

Geomagnetic index 

ATMPATH 

character 

“null” 

Path name for “atmosdat” atmospheric data file 

F10 

solar flux 
units 

0.0 

Daily 10.7-cmflux 

FlOb 

solar flux 
units 

0.0 

Mean 10.7-cm flux 

1ATMFL1 

integer 

1 

The GRAM initialization flag. 

For first vehicle number that uses GRAM 
=0, Do not initialize GRAM 
=1, Initialize GRAM 
Required for first call to GRAM 




For other vehicles using GRAM 
= 0, Do not update initial random number - will 
maintain same atmosphere density variability 
profile as another vehicle (correlated 
atmospheres) 

= 1 , Initialize atmosphere for new vehicle - no 
correlation between vehicle atmospheres 
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Input 

Symbol 

Type/ 

Units 

Stored 

Value 

Definition 

1ATMFL2 

integer 

0 

Random number seed - used when variable 
atmosphere mode (1ATMFL3=1) is used 

1ATMFL3 

integer 

0 

Perturbed atmosphere property options 
= 0, Mean density, temperature, pressure and 
speed of sound. 

= 1 , Perturbed density, mean temperature and 
pressure, and perturbed speed of 
sound 

1ATMFL4 

integer 

0 

Winds to use if NPC(6)=4 

= 0, nominal GRAM winds 

= 1 , nominal GRAM winds + perturbations 

IAUXPROFILE 

integer 

0 

Flag to indicate POST2 tables are used to define 
auxiliary profile. The auxiliary profile is used to 
override the default density, temperature, pressure and 
wind profiles. The use of this option does not allow 
for DUSTAU perturbations. 

= 0, Do not use POST2 tables for profile, use file 
define by PROFILE 

= 1, Use POST2 tables AUXTEMPT, AUXPREST, 
AUXDENST, AUXEWW1NDT, and AUXNSW1NDT 
to define auxiliary profile 

10PR 

integer 

0 

random output option 

1 = random output 

2 = none 

IDA 

integer 

0 

Day of the month 

NCEPHR 

integer 

5 

Code for UT hour of day if NCEP climatology is used: 
1=00 UT, 

2=06UT, 

3=12UT, 

4=1 8UT, 

5=all times of day combined, 

0 to use NCEP time-of-day based on input UTC hour 
(ihro) 

NCEPYR 

integer 

9008 

yly2 to use NCEP climatology for period-of-record 
(POR) from year yl through year y2 
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Input 

Symbol 

Type/ 

Units 

Stored 

Value 

Definition 

1HRO 

integer 

0 

Initial UTC (Greenwich) time hour 

1N1TPERT 

integer 

0 

Initial perturbations flag 
= 0 GRAM-derived random initial perturbations 
values 

= 1 User-selected initial perturbations 

ITHERM 

integer 

0 

Thermosphere model selection flag. 
= 1 MET (Jacchia) 

= 2 MS1S 
= 3 JB2006 

IUS 

integer 

0 

Unit number for atmosdat data 

IUG 

integer 

0 

Unit number for NCEP files 

1URRA 

integer 

0 

Unit number for Range Reference Atmosphere (RRA) 
data. 

= 0 none 

= xx actual unit number (recommend 42) 

1YR 

integer 

0 

4 digit or 2 digit year. If 2 digits are used 
lYR>56=19xx, 1YR<57 =20xx 

M1NO 

integer 

0 

Initial UTC (Greenwich) time minute 

MN 

integer 

0 

Month (1-12) 

PATCHY 

integer 

0 

Patchiness in perturbation model. 
= 0 no patchiness 
^ 0 patchiness 

PROFILE 

character 

“null” 

(Optional) auxiliary profile input file name. The 
auxiliary file is used to override the default density, 
temperature, pressure and wind profiles. 

RD1N1T 

percent 

0.0 

Initial density perturbation model (percent of mean). 
Used if 1NITPERT = 1 

RP1NIT 

percent 

0.0 

Initial pressure perturbation model (percent of mean). 
Used if INITPERT = 1 

RPSCALE 

decimal 

1.0 

Random density perturbation scale factor (0 = no 
perturbation) 0 < RPSCALE < 2. If the table 
RPSCALET is used, then the table value will override 
this input. 

RRAPATH 

character 

“null” 

Directory for RRA data. 
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Input 

Symbol 

Type/ 

Units 

Stored 

Value 

Definition 

RT1N1T 

percent 

0.0 

Initial temperature perturbation model (percent of 
mean). Used if 1N1TPERT = 1 

RUINIT 

m/s 

0.0 

Initial eastward wind velocity perturbation model. 
Used if 1N1TPERT = 1 

RVIN1T 

m/s 

0.0 

Initial northward wind velocity perturbation model. 
Used if 1N1TPERT = 1 

RWIN1T 

m/s 

0.0 

Initial upward wind velocity perturbation model. Used 
if IN1TPERT = 1 

SECO 

decimal 

0.0 

Initial UTC (Greenwich) time minute 

SITEL1M 

degree 

0.0 

Lat-Lon radius from RRA or PROFILE outside of 
which the RRA or PROFILE data are not used. 

S1TENEAR 

degree 

0.0 

Lat-Lon radius from RRA or PROFILE inside of 
which the RRA or PROFILE data are used with 1.0 
weighting factor. Between S1TENEAR and S1TELIM 
the weighting factor transitions from 1.0 to 0.0 
smoothly. 

S10 

decimal 

0.0 

Extreme Ultraviolet (EUV) index (26-34 nm) scaled to 
F10 units, i.e., 0.0 corresponds to S10=F10 

S10B 

decimal 

0.0 

EUV 81 -day center-averaged index scaled to FI 0B 
units (i.e., 0.0 corresponds to S10B=F10B 

XM10 

decimal 

0.0 

MG2 index scaled to F10 units (i.e., 0.0 corresponds to 
XM10=F10 

XM10B 

decimal 

0.0 

MG2 8 1 -day center-average index scaled to F 1 0 units 
(i.e., 0.0 corresponds to XM10B=F10) 


Table 5.5-2. POST2 Tables for the Earth-GRAM 2010 Model 


Input 

Type/ 

Stored 


Symbol 

Units 

Value 

Definition 

AUXDENST 

kg/m^ 

0.0 

Auxiliary profile of natural log of atmospheric density. 
Used if 1AUXPROF1LE = 1. 

AUXEWW1NDT 

m/s 

0.0 

Auxiliary profile of East/West wind. Used if 
1AUXPROF 1LE = 1. 

AUXNSW1NDT 

m/s 

0.0 

Auxiliary profile of North/South wind. Used if 
IAUXPROFILE = 1. 
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AUXPREST 

N/m^ 

0.0 

Auxiliary profile of natural log of atmospheric 
pressure. Used if IAUXPROF1LE = 1. 

AUXTEMPT 

degree K 

0.0 

Auxiliary profile atmospheric temperature. Used if 
1AUXPROF1EE = 1. 

RPSCALET 

decimal 

0.0 

Random density perturbation scale factor (0 = no 
perturbation) 0 < RPSCALET < 2. If the table 
RPSCALET is used, then the table value will override 
the RPSCALE input. 


Table 5.5-3. POST2 Outputs for the Earth-GRAM 2010 Model 


Output 

Symbol 

Type/ 

Units 

Definition 

ARMOLE 

percent 

Percentage of AR by volume in atmosphere. 

ATEM 

degree K 

Atmospheric temperature. 

AUXDENS 

kg/m^ 

Auxiliary profile density. Calculated if 1AUXPROFILE = 1. 

AUXEWWIND 

m/s 

Auxiliary profile East/West wind. Calculated if 
1AUXPROF ILE = 1. 

AUXNSWIND 

m/s 

Auxiliary profile North/South wind. Calculated if 
1AUXPROF1LE = 1. 

AUXPRES 

N/m^ 

Auxiliary profile pressure. Calculated if IAUXPROFILE = 1. 

AUXTEMP 

degree K 

Auxiliary profile temperature. Used if IAUXPROFILE = 1. 

CH4MOLE 

percent 

Percentage of CH4 by volume in atmosphere. 

COMOLE 

percent 

Percentage of CO by volume in atmosphere. 

C02MOLE 

percent 

Percentage of C02 by volume in atmosphere. 

CS 

m/s 

Speed of sound. 

DENS 

kg/nU 

Atmospheric density. 

DENS = DENS * GENTAB(DENKT) 

DENSMEAN 

kg/m^ 

Mean atmospheric density. 

DENSM3S 

kg/m^ 

3 sigma low atmospheric density. 

DENSP3S 

kg/m-'’ 

3 sigma high atmospheric density. 

DENSRAT 

decimal 

Ratio of density to mean density 

DENS76STAND 

kg/rrC 

1976 standard atmospheric density. 

EWWINDMEAN 

m/s 

Mean East/West wind velocity. Positive to the East. 

HEMOLE 

percent 

Percentage of HE by volume in atmosphere. 

HMOLE 

percent 

Percentage of H by volume in atmosphere. 

H20MOLE 

percent 

Percentage of H 2 0 by volume in atmosphere. 

MOLWE1GHT 

decimal 

Molecular weight of atmosphere 

NMOLE 

percent 

Percentage of N by volume in atmosphere. 

N2MOLE 

percent 

Percentage of N 2 by volume in atmosphere. 

N20MOLE 

percent 

Percentage of N 2 0 by volume in atmosphere. 

NSWINDMEAN 

m/s 

Mean North/South wind velocity. Positive to the North. 
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Output 

Symbol 

Type/ 

Units 

Definition 

OMOLE 

percent 

Percentage of O by volume in atmosphere. 

02MOLE 

percent 

Percentage of CE by volume in atmosphere. 

03MOLE 

percent 

Percentage of 0 3 by volume in atmosphere. 

PRES 

N/m^ 

Atmospheric pressure. 

PRESMEAN 

N/m^ 

Mean atmospheric pressure. 

PRES76STAND 

N/m^ 

1976 standard atmospheric pressure. 

TEMPMEAN 

degree K 

Mean atmospheric temperature. 

TEMP76STAND 

degree K 

1976 standard atmospheric temperature. 

VERTW1NDMEAN 

m/s 

Mean vertical wind velocity. Positive down. 

W1NDEW 

m/s 

East/West wind velocity. Positive East. Used if NPC(6) =4 

W1NDNS 

m/s 

North/South wind velocity. Positive North. Used if NPC(6) =4 

W1NDV 

m/s 

Vertical wind velocity. Positive down. Used if NPC(6) =4 

ZHPRES 

km 

Atmospheric pressure scale height. 

ZHRHO 

km 

Atmospheric density scale height. 


5.5.2 Earth Gravity Model 96 

The EGM96 model was formatted for use with the POST2 gravity model. The EGM96 
(complete to degree and order 360) was generated by NASA and other Government agencies. 
This model is being used with various NASA simulations evaluating currently proposed Earth 
entry systems. Various data sets including satellite tracking and altimetry were used to generate 
this model. Standard POST2 inputs and outputs for gravity models are used; the spherical 
harmonic gravity model is invoked in POST2 using NPC(16)=7, and the file containing the 
sectoral and tesseral data is identified in the input (including system directory path) using the 
POST2 variable GRAVDATA. A test input file was generated using EGM96. 

5.6 Attitude Control Models 

In an effort to balance faster executing and more easily developed 3 DoF simulations (versus the 
slower, higher fidelity 6 DoF simulations), control models that address vehicle attitude change 
are used. A pseudo controller model was included in Phase 1 which used a method of emulating 
a 6 DoF attitude control system for bank angle modulation. Additionally in Phase 1, the natural 
aerodynamic trim points in pitch and yaw to determine the angle of attack (AO A) and sideslip 
angle that the vehicle orients to at any time during the atmospheric entry were included to 
provide 3-axis vehicle control. For Phase 2, another method of 3-axis vehicle control based on 
the bank angle pseudo-controller was developed. Also, a utility function to provide a second 
order actuator model for use with models within POST2 (e.g., engine gimbal or sensor slewing 
motion) is included. Both of these models were incorporated in the simulation and their 
implementations are described in the following subsections. 
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5.6.1 Pseudo-Controller for 3-axis Attitude Control 

The program contains a pseudo-controller module that emulates the performance of a feedback 
control system. The three-axis pseudo-controller (3APC is a generalization of the bank angle 
pseudo-controller (BPC) and allows a 3 DoF simulation to model some of the dynamic attitude 
behavior of a 6 DoF trajectory. The 3APC can be used for 3-axis attitude control of any 3 DoF 
vehicle trajectory. However, if the vehicle being modeled uses bank angle to modulate lift, then 
the BPC will be a more straightforward implementation. Note that this module does not solve 
the rotational equations of motion and does not require moments of inertia. Also, the 3APC 
module treats each axis independently and should not be used for systems with high angular rates 
or other significant cross-axis coupling. 


The 3APC expects three commanded angles from the guidance system: roll, pitch and yaw. The 
3APC attempts to achieve these commanded angles by applying angular accelerations under the 
constraint of second-order dynamics. Maximum allowed values for the angular acceleration and 
angular rate are given as inputs: 


<t>- 


+ 


(j)d t+ <j)d ? 


The 3APC module maintains internal states of angle and angle rate for each axis, which are 
integrated within the model, not using the POST2 generalized integration procedure. Angular 
rate and acceleration can be adjusted based on current residual angle error (i.e., difference 
between commanded and actual angle) and when the actual angle overshoots the command. 

Several options are available to model behavior of the 3APC when the vehicle overshoots the 
command. Since the 3APC bases the angular acceleration, it uses on the current angle and the 
commanded angle, overshoots should only occur due to a change in the command. The 
“normal” option is to use all of the vehicles deceleration capability to reach the commanded 
angle. The “no overshoot” option allows the vehicle to exceed the maximum angular 
acceleration when the vehicle passes the commanded angle in order to prevent an overshoot. 

The “no wrong way” option ensures that the angle only moves toward the commanded angle. If 
the command changes in such a way that the current angular rate moves the vehicle away from 
the current command, the “no wrong way” option forces the angle to the new command. Finally, 
the “perfect” controller option instantly sets the angle to the command and the angular rate to 
zero. 

The angle command and the angle direction command are provided from the guidance or the 
input deck. There are several options for angle direction and for how the controller handles 
overshoots. The angle direction flag can be set to: under, left, shortest, right, or over. “Under” 
means the vehicle will rotate toward an angle of 180 degrees. “Left” means the vehicle will 
rotate to the left. “Shortest” means bank in the direction that minimizes the absolute angle error. 
“Right” means the vehicle will rotate to the right, and “over” means through an angle of 0 
degrees. 
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No limitations are placed on the size of the maximum acceleration or rate that are input. If large 
enough values are used, then the 3APC will operate as a near instantaneous angle controller. 
Once input, these limits are used to limit the maximum values in an absolute value sense (i.e., + 
maximum value are the limit boundaries used). 

Inputs and outputs for the POST2 implementation of the 3APC are provided in Tables 5.6-1 and 
5.6-2, respectively. To use the 3APC, the relative euler angles steering option must be used 
(IGUID(1)=2), and the angle polynomial with constant term from input (IGUID(4) =1 or 
IGUID(9) =3*1, depending on value of IGUID(2). 


Table 5.6-1. POST2 Inputs for the 3APC Model 


Input 

Symbol 

Units 

Stored 

Value 

Definition 

CTRL_ROL_CMD , 
CTRL_PIT_CMD , 
CTRL_ Y A W_CMD 

degree 

0.0 

Commanded angle about the X (ROL), Y (PIT), 
and Z (YAW) body axes. 

CTRL_ROL_DT, 
CTRL_PIT_DT, 
CTRL_ Y A W_DT 

seconds 

0.0 

Pseudo-controller update cycle time for the X 
(ROL), Y (PIT), and Z (YAW) body axis 
controller. 

CTRL_ROL_IDIR_FLAG, 
CTRL_PIT_IDIR_FLAG, 
CTRL_Y A W_IDIR_FLAG 

integer 

0 

Flag to control maneuver direction for the X 
(ROL), Y (PIT), and Z (YAW) body axis 
controller 

= -2, bank through 1 80 degree (underneath) 
= - 1 , bank left 
= 0, go shortest distance 
= 1 , bank right 

= 2, bank through 0 degree (over) 

CTRL_ROL_ 

IPSEUDO_FLAG, 

CTRL_PIT_ 

IPSEUDO_FLAG, 

CTRL_YAW_ 

IPSEUDO_FLAG 

integer 

0 

The Pseudo-controller selection flag for the X 
(ROL), Y (PIT), and Z (YAW) body axis 
controller. 

= 1 , Use Pseudo-Controller 

CTRL_ROL_MAXACCEL, 
CTRL_PIT_MAXACCEL, 
CTRL_Y A W JVIAXACCEL 

deg/s 2 

5.0 

Maximum angular acceleration used by the 
pseudo-controller flag for the X (ROL), Y (PIT), 
and Z (YAW) body axis 

CTRL_ROL_MAXRATE, 
CTRL_PIT_MAXRATE, 
CTRL_Y A W_M AXRATE 

deg/s 

20.0 

Maximum angular rate used by the pseudo- 
controller flag for the X (ROL), Y (PIT), and Z 
(YAW) body axis 
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Input 

Symbol 

Units 

Stored 

Value 

Definition 

CTRL_ROL_OSMODE, 

CTRL_PIT_OSMODE, 

CTRL_YAW_OSMODE 

integer 

0 

Overshoot mode selection flag for the X (ROL), Y 
(PIT), and Z (YAW) body axis 

= 0, normal 

= 1 , no overshoot 

- 2, no wrong way 

= - 1 , perfect 

IGUID(l) 

integer 

0 

Type of steering (guidance) desired. 

=0, AOA, sideslip, and bank. Also 
input values for IGUID(2) and IGUID(3) or 
IGUID(6), IGUID(7), and IGUID(8). 

IGUID(3) 

integer 

0 

A flag to specify the steering option when 
commanding all channels simultaneously using 
aerodynamic AOA, sideslip, and bank angle. Must 
use option 1 with BPC. 

= 1 , Command AOA, sideslip, and bank as third 
order polynomials with the values of the constant 
terms of the polynomials are the input values. 

IGUID(8) 

integer 

0 

Steering option flag when using separate channel 
for bank angle. Must use option 1 with BPC. 

= 1 , Command bank angle as third order 
polynomials except that the constant terms of the 
polynomials are the input values. 


Table 5.6-2. POST2 Outputs for the 3APC Model 


Output 

Symbol 

Type/ 

Units 

Definition 

CTRL_ROL_CMD , 
CTRL_PIT_CMD , 
CTRL_Y A W_CMD 

degree 

Angle command for the X (ROL), Y (PIT), and Z (YAW) body 
axis 

CTRL_ROL_DOT, 
CTRL_PIT_DOT, 
CTRL_Y AW_DOT 

deg/s 

Angular rate about the X (ROL), Y (PIT), and Z (YAW) body 
axis 

CTRL_ROL_PC 1 , 
CTRL_PIT_PC 1 , 
CTRL_Y A W_PC 1 

degree 

X (ROL), Y (PIT), and Z (YAW) body axis angle polynomial first 
coefficient - generated by pseudo-controller 
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5.6.2 Second Order Actuator Model 

A general purpose second order actuator model has been developed for use internally with 
POST2 models (e.g., engine gimbal or sensor slewing). The routine, actu8, models second order 
dynamic system of the form: 



+ 2 4® 


dx 

dt 


+ co 2 x = u 


(5.6-1) 


where x is the displacement of the system; Z, is the damping; cd is the natural frequency; and u is 
the input to the system. The model expects to be called at every time step and integrates the 
1 internal state variables across a time interval that must be supplied. The user provides the 
actuator damping and natural frequency, the update cycle time, the position and rate limits, and 
the current position command. 

Inputs to the routine are: 

• System initial position 

• System initial rate 

• Time step 

• Actuator damping ratio 

• Actuator natural frequency 

• Commanded position 

• Plus and minus rate limits 

• Plus and minus displacement limits 

Outputs from the system are: 

• Final position after time step 

• Final rate after time step 

Note that the final position and rate from one call are used as the initial position and rate for the 
next call. This feature allows the same actuator model code to be used for multiple actuators, at 
the expense of external storage for the state variables. 

6.0 Observation and NESC Recommendation 

The following observation and NESC recommendation relate to technical aspects of this 
assessment: 

0-1. The validated and documented subsystem models, test cases, and scripts developed 

during Phases 1 and 2 of this assessment, and added to the POST2 simulation framework, 
have increased the Agency’s ability to rapidly evaluate EDL characteristics in systems 
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analysis studies, preliminary design, mission development and execution, and time- 
critical assessments. 

R-l. The custodians of the POST2 simulation framework, the Atmospheric Flight and Entry 
Systems Branch at LaRC, should retain and maintain the REDLAS models, test cases, 
and scripts for use by current and future programs. ( 0 - 1 ) 

7.0 Other Deliverables 

There are no other deliverables. 


8.0 Definition of Terms 

Corrective Actions Changes to design processes, work instructions, workmanship practices, 
training, inspections, tests, procedures, specifications, drawings, tools, 
equipment, facilities, resources, or material that result in preventing, 
minimizing, or limiting the potential for recurrence of a problem. 

Finding A conclusion based on facts established by the investigating authority. 

Lessons Learned Knowledge or understanding gained by experience. The experience may 
be positive, as in a successful test or mission, or negative, as in a mishap 
or failure. A lesson must be significant in that it has real or assumed 
impact on operations; valid in that it is factually and technically correct; 
and applicable in that it identifies a specific design, process, or decision 
that reduces or limits the potential for failures and mishaps, or reinforces a 
positive result. 


Observation A factor, event, or circumstance identified during the assessment that did 

not contribute to the problem, but if left uncorrected has the potential to 
cause a mishap, injury, or increase the severity should a mishap occur. 
Alternatively, an observation could be a positive acknowledgement of a 
Center/Program/Project/Organization’s operational structure, tools, and/or 
support provided. 


Problem 


The subject of the independent technical assessment. 


Proximate Cause The event(s) that occurred, including any condition(s) that existed 
immediately before the undesired outcome, directly resulted in its 
occurrence and, if eliminated or modified, would have prevented the 
undesired outcome. 
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Recommendation An action identified by the NESC to correct a root cause or deficiency 

identified during the investigation. The recommendations may be used by 
the responsible Center/Program/Project/Organization in the preparation of 
a corrective action plan. 


Root Cause 

One of multiple factors (events, conditions, or organizational factors) that 
contributed to or created the proximate cause and subsequent undesired 
outcome and, if eliminated or modified, would have prevented the 
undesired outcome. Typically, multiple root causes contribute to an 
undesired outcome. 

9.0 Acronyms List 

3APC 

Three Axis Pseudo-Controller 

AMA 

Analytical Mechanics Associates, Inc. 

AOA 

Angle-Of-Attack 

ARC 

Ames Research Center 

ATK 

Alliant Techsystems, Inc. 

BPC 

Bank Angle Pseudo-Controller 

CFD 

Computational Fluid Dynamics 

eg 

Center of Gravity 

DFRC 

Dryden Flight Research Center 

DoF 

Degree of Freedom 

EDF 

Entry, Descent, and Fanding 

EDF-SA 

EDF Systems Analysis 

EGM 

Earth Gravity Model 

EKE 

Extended Kalman Filter 

EUV 

Extreme Ultraviolet 

GN&C 

Guidance, Navigation, and Control 

GRAM 

Global Reference Atmosphere Model 

GSFC 

Goddard Space Flight Center 

IMU 

Inertial Measurement Unit 

JSC 

Johnson Space Center 

KSC 

Kennedy Space Center 

LaRC 

Langley Research Center 

LOX 

Liquid Oxygen 

mEKF 

multimode EKF 

MER 

Mars Exploration Rover 

MSFC 

Marshall Space Flight Center 

MSL 

Mars Science Laboratory 
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n 2 

NCEP 

nd 

NESC 

NIA 

NRB 

0 2 

0 3 

POR 

POST2 

REDLAS 

RRA 

TPC 


nitrogen 

National Centers for Environmental Prediction 
non-dimensional 

NASA Engineering and Safety Center 
National Institute for Aerospace 
NESC Review Board 
oxygen 

the molecular formula for ozone 
period-of-record 

Program to Optimize Simulated Trajectories II 
Rapid EDL Analysis Simulation 
Range Reference Atmosphere 
Terminal Point Control 
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Appendix A. Multi-Mode Extended Kalman Filter 

Work in this Appendix was performed under contract NNL07AA00B, 

Task Order NNL07AM00T 
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1 Introduction 

The main objectives of the project are: 

- Develop a multimode Kalman navigation filter (mKNF) architecture for application to entry, de- 
scent, and landing (EDL) using a standard set of measurements including altimeters, velocimeters, 
star trackers, and inertial measurement units. The main computational unit of the navigation archi- 
tecture is die extended Kalman filter algorithm employing a dual-inertial state implementation with 
integrated position, velocity, and attitude. The filter will be readily re-configurable to accommodate 
differing mission scenarios (such as lunar entry versus Mars entry). 

- Develop sensor models as required to support the mKNF and requisite mathematical specifications 
documentation. The error models will include systematic and random errors applicable to the in- 
dividual sensors, including biases, misalignments, scaling factors, and random noise. 

- Develop and execute test scripts to verify and validate die mKNF. This includes specifying test proce- 
dures, creating run matrices, reviewing the simulation outputs, creating run summaries, and com- 
menting on the results. 

To accomplish the tasks, algorithms were packaged in compact units of code with input/ ouput tightly 
controlled using data structures. Each individual unit of code would be entirely re-entrant; hence the 
use of global variables and persistent variables is disallowed. Two main hierarchical structures were de- 
veloped: the Multi-mode EKF Hierarchical Structure (MHS) and the Sensor Hierarchical Structure (SHS). 
The MHS shown in Figure 1 enables the mEKF to be a reusable unit from multiple calling programs. All 
of die inputs and outputs to die mEKF shown in Figure 1 are data structures with definitions described in 
Section 4. 


i 1 



Figure 1 : Multe-mode EKF Hierarchical Structure 

Similarly, a sensor hierarchical structure (SHS), shown in Figure 2, is proposed to represent a convenient 
vehicle for representing the sensors required for EDL. The SHS is comprised of four major substructures: 
Universe, Nucleus, Identity, and Errors. 
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Figure 2: Sensor Hierarchical Structure 


The Universe structure provides the parameters required by the filter that describe the external en- 
vironment. The Nucleus structure describes the key elements of the sensor around which the model is 
constructed including, for example, the rate (internal), output rate, on/ off events, and special process- 
ing (hardware specific) parameters. Currently, the Nucleus structure models the velocimeter, spherical 
altimeter, attitude camera, and IMU. The Identity structure describes the physical aspects of the sensor 
such as, for example, the spacecraft sensor layout (sensor locations, sensor reference frame orientations, 
etc.), and beam locations (for the velocimeter and altimeter). The Errors structure includes the random 
and systematic errors associated with each sensor, such as biases, noise strengths, misalignments, and 
time constants. Details of the SHS are provided in Section 3 for the IMU and in Section 5 for the altimeter, 
velocimeter, and attitude camera. 
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2 Extended Kalman Filtering Basics 

2.1 Architecture 

An extended Kalman filter (EKF) is a model-based algorithm. Hie ability of the EKF to provide accurate 
estimates of the state of the system depends explicitly on the quality of the model. In the development of 
the EKF we assume a nonlinear system model is assumed to be of die form 

x(i) = f(x(t),*) + M(t)w(t), (1) 

where x(t) c R n is the state of the system, f(x(t),i) c is a sufficiently differentiable nonlinear system 
model, and w(t) € R p is the process noise, which is assumed to have die properties 

E {w (l)} = 0 and E {w(t)w T (r) } = Q ap€ cS(t - r), 

where E { } is the expectation operator and 5(t r) is the Dirac delta iunction. The process noise is 
assumed to be a zero-mean, white-noise process with a constant-power spectral density given by Q spec e 
R pxp , which is an input to the EKF. Furdiemiore, M(f.) £ R nxp is the process noise mapping matrix. 

The nonlinear measurement model is assumed to be of the form 

y k = bk{x k ) + v k (2) 

where die subscript "k" denotes a discrete lime measurement at time t ^ Furdiermore, h k(x-k) € R m is 

a sufficiently differentiable nonlinear measurement model evaluated at the state x* = x(ffe), and v*. € R m 
is die measurement noise, assumed to be a white noise sequence with 

E{v*} = 0 and E{v k vl,} =R k 6 kk ,, 

where k and k' me two different discrete measurement times and 6kk' is tfie Kronecker delta. Additionally, 
it is assumed that the process noise and the measurement noise are not correlated in time, or, that 

E{w(f.)v£} =0 Vt,t k . 

Letting x(i) denote die true state at some time l and dial die estimate of the true state is denoted by 
x(t), define the estimation error as 

e(t)=x(t)-x(t). (3) 

By design, die EKF algoridim creates an unbiased state estimadon, so dial 

E(e(t)} = 0 V t. (4) 

This requires that the initial state estimate x(*o) is such that the initial estimation error is zero mean. This 
is accomplished by setting x(i 0 ) = E {x(£ 0 )}. We also define die state estimation error covariance as 

E{e(t)e T (t)} = P(t) Vt. (5) 

The EKF is a recursive algorithm comprised of two main phases: propagation and update. Figure 3 
illustrates die propagation phase and die update phase. The state estimate and the state estimation er- 
ror covariance are propagated forward to the time a discrete measurement becomes available at time t k . 
At this time, the measurement is used to update botii die state estimate and the state estimation error 
covariance. The values immediately preceding an update are called as a priori and denoted widi a super- 
script while the values immediately following an update are called as a posteriori and denoted widi a 
superscript The state estimation error covariance is a measure of uncertainly in the state estimate, so 
w’e expect it to be reduced (on average) as measurements are processed during die update phase. In die 
mF.KF, the external sensors (cunently) include the attitude camera, spherical altimeter, and velocimeter. 
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Altimeter, Veloci meter, and Attitude Camera 



Figure 3: EKF Timeline. 


It is during the update phase that the model of the sensor plays a key role. During the propagation phase 
there are no external measurements, so we expect (on average) that the state estimation error covariance 
increases (in other words, our uncertainty in the state estimate increases). During the propagation phase, 
the model of the dynamics of the system plays a key role. The inertial measurement unit (I MU) assists in 
the propagation of the state estimate by measuring and accounting for the nongravitational accelerations 
acting on the system. 


2.2 Propagation 

Taking the expected value of Eq. (1) yields 

= E {*(*)} = E{f(x(*),i)} + M(i)E{w(i)}. 

Expanding f(x(<), t) in a Taylor series expansion about x(t) and neglecting higher-order terms yields 

x{t) = E{f(x(t),t) + F(x(t),t) (x(t) - x(t))} + M(t)E (w(t)} , (6) 


where 




<9f(x(i),i) I 
^*(0 lx(t)=*(t) 


Eq. (6) can be simplified as 

x(t) = f(x(t),i) + F(x(t),i)E{(x(i) - x(£))} + M(<)E{w(<)} , 


Noting that the EKF is an unbiased estimator (E {(x(i) - x(£))} = 0) and that the process noise is assumed 
to be zero-mean (E (w(£)} = 0), the differential equation for the estimated state can be written as 

i(t) = f (*(*),*). (7) 

Using die estimation error in Eq. (3), and taking die time derivative yields 

e(t) =x(t)~ x(t). (8) 

Substituting x(£) from Eq. (1) and x(t) from Eq. (7) into Eq. (8) it follows tiiat 

e(t) = f (x(t),t) + M(t)w(t) - f(x(t),t). 
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Using a Taylor series expansion of the true dynamics, represented by f(x({), i), about the estimated state 
of the system, represented by f (x(t), t), and neglecting the higher-order terms yields 

e(t) = F(x(i),t) (x(i) x(t)) + M(t)w(t). 

Tlterefore, making use of the estimation error as defined in Eq. (3) and, after algebraic manipulation, we 
obtain die estimation error differential equation 


e(t) = F (x(t), f)e(i) + M(<)w(f ). (9) 

the solution of Eq. (9) is 

e(f) = $(i,fo)e(J 0 )+ / $(i,T)M(7-)w(r)£i7-, 

J to 

for all t > to, where the state transition matrix <£(#, to) corresponds to F(x (£)» £)> and satisfies die matrix 
differential equation 

$(Mo)=F(x(0,*)*(Mo), 

subject to die initial condition 

= Inxn- 

Frora the definition of the estimation error covariance in Eq, (4), it follows that 


P(f) = *(t,fo)P(«o)<b T (f, to) + /* $(f,r)M(r)Q, peP (7-)M :r (r)$ 7 '(i, r)dr. 

Jto 

Define the matrix Q(t) e R nxn as 

Q (<)= / *(«, r)M(r)Q, 1> «(r)M T (r)® T (t, r)dr. 

JtQ 

Taking the time derivative of Q(£) yields 

Q(i) = F(x(i),i)Q(f) + Q(t)F r (x(t),t) + M(0Q Jf «c(i)M T (f), 

with the initial condition 

Q(*o) = Onxn- 

To propagate the state estimate and the state estimation error covariance on the interval tk~i <t<tk, 
we numerically integrate 

= F (x(t),t)*(t,i*-i) 

Q(() = F(x(<),t)Q(f) + Q(f)F T (x(t),f) + M(/)Q 3p «(i)M T (t) 

*(*) = f (*(*)•*> 

from t - tk-i to t, — t-k with initial conditions 


h-l) = I„xn 

Q(ffc-i) = 0„x» 

X(t*_!) = X+_! (10) 

to obtain Q(£*) and i). After integration, die state estimation error covariance is mapped forward 

via 

P * - **PT-1*» + Q* 

x* (11) 

where Q k = Q(h-) and <!>,, = «(t*, d*~i). 
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2.3 Update 

Taking the expected value of both sides of Eq. (2) yields 

y * =E{y fc } =E{h fc (x*)} + E{v fc }. 

Using a Taylor series expansion of hj(xfc) about the a priori state estimate, x s: , and neglecting higher- 
order terms it is seen that 


y* = E{h*(x t ) +H 4 (x t ) (x* -* t )} + E{v*}. 

where 1 1, (x k ) is the measurement sensitivity matrix, defined to be 


«*(**) 


dh t (x t ) | 
dx k | 


lx*=it t 


Eq. (12) can be written as 


( 12 ) 


h = h k (x k ) + H k (x k )E{(x k -x k )} + E{v k }. 

Since tlie EKF is an unbiased estimator (E { (x* - x k ) } = 0) and the measurement noise is assumed to 
be zero-mean ( E {v*} = 0), the estimated measurement is found to be 


y k = h*(x t ). 


Note that both the measurement and the estimated measurement are evaluated just prior to the mea- 
surement update. Defining the measurement residual, r*, and the Kalman gain matrix, K*., 10 be 


= yfc — y k (13) 

K/e = P- k Hl(Sr k )W- k \ (14) 

where W k is the measurement residual covariance matrix, 

W* = H*(**)P*Hi[(**) + R*, (15) 

the state estimate and the state estimation error covariance update are then calculated as 

** = + K k r k (16) 

Pt = [l„x„ — KfcHjt(xJ)] [l nxn — KfcHjt(xJ)] r + KfcRfcKfc. (17) 


In the EKF development, we must specify the dynamics model f(x(t), t), the measurement model, hfc(xfc), 
the process noise matrix, Q spec , the process noise mapping matrix, M(<), and die measurement noise co- 
variance Rfc. To execute the EKF we need to specify the initial conditions x(£ 0 ) = E{x(f 0 )} and P 0 — 
E | (x(£ 0 ) - x(£o)) (x(£ 0 ) - x(t 0 )) T | • As die measurements become available at t k , we can propagate die 
state estimate and state estimation error covariance, xf_j and Pj_ lf , respectively, from t k _i to t k to ob- 
tain x* and P^ , and use die measurement information to update die state estimate and state estimation 
error covariance to obtain xj and PJ. This process continues as the measurements continue to become 
available. Overall , die EKF architecture is algorithmically straightforward. The challenge is to determine 
a viable set of states and associated models for the dynamics and sensors so diat die EKF can be imple- 
mented in real-time. In this work, we also strive to create a structured architecture that permits the EKF 
to be readily reconfigured for different sensors and scenarios. 
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3 Mathematical Model of the Flight Dynamics 

3.1 Dynamical Model in Continuous-time 

The dynamics of the spacecraft during EDL occur naturally in continuous-time. However, since the 
EKF executes on a digital computer, the dynamics modeled in the estimation algorithm occur naturally 
in discrete-time. Also, the I MU and external sensors provide measurements at discrete times. The strat- 
egy here is to present the dynamics of the spacecraft during EDL in continuous- time and then with ap- 
propriate assumptions create a discrete-time model equivalent. One important aspect of the approach 
followed here is to utilize the IMU to provide a measurement of all the non-gravitational accelerations 
(thrust, aerodynamics, etc.) acting on the spacecraft rather than attempting to model their associated 
dynamics explicitly. 

In the mEKF we navigate the IMU (not the center of gravity). The total acceleration is comprised of 
accelerations due to gravity and accelerations due to non-conservative forces, such as thrust and aerody- 
namics. The acceleration due to gravity is a function of the position of the center-of-gravity of the vehicle, 
where 

r c 3 = r L.Jt) + T b(4i) r c s /imu- 

The IMU is not generally co-located with the CG of the vehicle; instead, it is offset at an offset position 
from the CG, denoted by rj 9 / imu ,v^ ere the superscript b denotes body reference frame. The transforma- 
tion matrix from body reference frame to inertial reference frame, denoted by T^(q^), is given by 

TH«?)= [(9 2 + q T q)l3*3-29[qx]+2[qxf] T €R 3 * 3 , (18) 

where 



is the attitude quaternion representing the inertial reference frame relative to the case reference frame. 
The continuous equations of motion are given by 


HmuW = v imu ( 1 ) 

+ Ti(q^T b c a c ng (t) 

= JO}/, (*)®«f(9 

where Tj is a known constant matrix and u> is the pure quaternion given by 


&b/i = 



(19) 

( 20 ) 

( 21 ) 


The non-conservative forces (accelerations denoted by a c ng {t) in Eq. (20)) are sensed by the IMU. The IMU 
then becomes an instrumental element of die state propagation since it essentially replaces a mathemat- 
ical model of the non-conservative forces. 

Note that the quaternion of rotation is subject to a unity norm constraint or 

g = |4l = VV 2 + q T q = l. 

Quaternion multiplication (i.e. the composition of rotations) for any two arbitrary quaternions, p and f 
is given by 

p<8> f = | ^ J ® 

10 


pr + pr - p x r 
pr — p T r 


( 22 ) 
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where the symbol <S> is exclusively used to denote quaternion multiplication. In the special case that the 
quaternion is used to represent a small rotation, denoted by Sa, the small angle quaternion, denoted by 
Sq, is given by 



\Sa 1 

1 J ' 


(23) 


Substituting Eq. (23) into Eq. (18) yields (to first order in <5a), the rotation matrix, <5T, representing the 
small angle rotation 

<5T = I 3x 3 — [Sax ] , (24) 


where 


[<k*x] 


0 — 6 a$ 6 a ,2 

60,3 0 — 6 a\ 

— 8 o 2 0 


Therefore, given a set of small angles which may represent a deviation in attitude, the corresponding 
quaternion is given via Eq. (23) and the corresponding rotation matrix is given via Eq. (24). This is impor- 
tant because the attitude representation using quaternions increases the complexity of the estimation 
process since the attitude update occurs through multiplication rather than through addition. The mEKF 
will utilize a multiplicative update for the attitude estimation based on small angle approximations. 


3.2 Inertial Measurement Units 

The strapdown I MU accelerometers and gyros produce measurements corrupted by random and sys- 
tematic errors. The accelerometer direcdy measures the non-gravitational accelerations acting on the 
spacecraft. We assume that the accelerometer package outputs a measure of the spacecraft change in 
velocity due to non-gravitational accelerations (expressed in the IMU case frame) over the time interval 
tk~ i to tk denoted by Av^ k . We model the computation of the true Av^ k via 

rt k 

^true.k = dr. 

Jtk-l 

The actual measurement of acceleration is corrupted by systematic and random errors 
errors we assume a model of the form 

J&(t) = (I sx3 + r a )(I 3 x3 + S„) (a c (t) + 0 a + i t ) , (26) 

where 



0 

7a« 

7a* y 



0 

0 

r„ := 

-7... 

0 

7 « v * 

, S a := 

0 

S “v 

0 


7 o sv 

—7 a« 

0 


0 

0 

s a t 


and 7 „ = ('fa,, , 7<.,„>7«„>7«„,7a„,> 7a,, ) T are nonorthogonality and axes misalignment errors, /9„ e R 3 
is the bias in the accelerometer, s a = (s 0l , s„ v , s a , ) T are scale factor errors, and e R 3 is a zero-mean, 
white noise process with E [£, = 0 and 

E J — E 6 (t — t) where S(t — r) — 0 when t ^ r and S(t — r) = 1 when t — r. 

The nonorthogonality and axes misalignment errors, scale factor errors, and bias parameters are all mod- 
eled as zero-mean random constants. We characterize the random constants with E [-yj = 0 ,E [sj = 0, 
E \p a ] = 0 and the given matrices E [ 7 ^ 7 ^] . E [s„s^] , and E ] • If we assume that the various errors 
are small, then to first-order 


(25) 

. To model the 


(Isx3 +ra)(l 3 x 3 + S a ) ft; I 3 X 3 + T a + S Q . 
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Defining 

a u r„ + s„ 

yields the accelerometer measurement model 


as, (0 = (Isxs + A„) (a'(t) + j8 a + € t ) , 


(27) 


(28) 


Typically, the IMU error parameters (such as the bias covariance and noise spectral density) associated 
with the continuous- time acceleration error model in Eq. (28) are provided as part of the hardware specifi- 
cations. The challenge here is to develop an error model that reflects the output of the IMU given in terms 
of integrated acceleration, understanding dtat there may be internal computations occurring within the 
IMU unit and that the IMU outputs are available only at discrete times. To that end, model the measured 
Avjj, k at time t* via 

Av5,jk= r a^(r)dr (29) 

■>tk - 1 

where a ((,(<) is given by Eq. (28). The associated accelerometer measurement model, AvJ,*, is repre- 
sented by 

A< 1>t = (I + A 0 )(A + b„ + i k ). (30) 

where A a (~/ a , s 0 ) are noriorthogonaliiv and axes misalignment errors and scale factor errors given as 
above in Eq. (28), b„ € E 3 is the bias in the accelerometer and is a function of the integration time, and 

as 


£ k £ K 3 is a zero-mean, white noise sequence. From Eq. (30), compute A v c mk < 

Av m,k = Av ?r ue ,fc + D ( Av m,*) S « + N ( Av m,t)7« + K + t k 
where we have utilized the approximation (I + A a ) -1 ~ I - A„ and 


f-K Av m,fc) ’ 


where Avf, 


(31) 


r a v* 

0 

0 


- Av z 

Aljy 

0 

0 

0 

0 

! 0 

Av y 

0 

. N(Av^) := 

0 

o' 

Av z 

— Ati-g 

0 

0 

0 

Cl 

A «z 


0 

0 

0 

0 

- A »» 

Av x 


v m j, = ( A 'v : A v y ,Avg) . In the current version of the mlKE we estimate only tfie IMU ac- 
celerometer bias, b 0 , and treat s„, 7 o, and £ k as contributing to the process noise of the dynamics. Hence, 
we have the model 

Av m,k = A < rue ,fc + bo + w ai ,, fc (32) 

where wa,,^ - D(Av‘ n i )s„ + N(Av' n i ) 7 „ + ( k and we assume that E - W ^„6 kj where 

6k j ~ 0 when k ^ j and ojy = 1 when k = j and we choose Wa„ appropriately. 

The situation with die gyros is similar to that of the accelerometers. The gyro measures the angular 
velocity of the spacecraft; however, the gyro package outputs a measure of the spacecraft attitude change 
over the time interval tk - 1 to tk denoted by A 0 c mk - The measurement of the angular velocity vector is 
corrupted by random biases, errors due to scale factor uncertainties, errors due to nonorthogonality and 
axes misalignments, and random noise. The gyro error model can be formulated as 


A 0?™, k = j dr 

Jt k - 1 

and the corresponding measurement model is 


(33) 


a c ,«- r 

jtk-i 


dr 


from which it follows that 


A ®m,k — (13x3 + r 9 )(I 3 x3 + S 9 )(A0J,. u<!ii + b g + f) k ) , 


(34) 
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where b p is hie gyro bias, S g is the gyro scale factor matrix, T g is die gyro nonordiogonality and axes 
misalignment matrix, rj k is a white sequence, and where 



s s. 

0 

0 


0 

7 s„ 

-7 g„ 

9 : — 

0 

% 

0 

. r 3 := 

“7*. 

0 

7w. 


0 

0 

%. 


. ~1g.y 

~l9.r 

0 


and 7 fl = ( 7 ^ , 7 9xv , 7 9zi , y gxz , y g . cw , y 9vx ) T are nonordiogonality and axes misalignment errors, h g € R 3 is 
the bias in die gyro, s g = (s 9x , s 9v , sg,) 1 are scale factor errors, and tj k € R 3 is a zero-mean, white noise 
sequence. The nonordiogonality and axes misalignment errors, scale factor errors, and bias parameters 
are all modeled as zero-mean random constants. To first-order, we have 


(13x3 + S 5 )(I 3 *3 + T g ) I 3n ;3 -f + T s . 

Hence, Eq. (34) can be written in the form 

&Qm,k — (^3X3 + &g)(&&true,k + + Wk)’ (3®) 

where A g :—S g + T g . From Eq. (36), compute 

A 0 c nitk = A r V(AO c m>k )s g 1 N(A0 c m>k )y g 4 b g 4 *?*, (37) 

where we liave utilized the approximation (I + A & ) _1 ~ I — A g and 



' A0 X 

0 

0 1 

<3* 

<1 

1 

A 0,j 

0 

0 

0 

0 ‘ 


0 

AS,, 

0 | , N(A0^ iJfc ) := 

0 

0 

A6 r , 

-a e x 

0 

0 


0 

0 

A t) z J 

0 

0 

0 

0 

-A0 y 

<1 


where A 0 c m k — (A 9 Xi A 0 y , A 0 Z ) . in die current version of die mEKE, we estimate only die 1MU gyro 
bias, b 9 , and treat s g , y g , and rj k as contributing to the process noise of the dynamics. Hence, we have the 
model 

A 0 r mk — A 0 c true _ k + b s + w^,k, (38) 

where = D(A0' lj; .)s 9 + N(Ae^ !>i ) 7 9 + T) k , and we assume that E = W ae 4j where 

Sfcj = 0 when k ^ j and 6 k j = 1 when k = j and we choose Waj appropriately. 

3.2.1 Accumulated Integrated Accelerometers and Gyros 

Consider die scenario in which an IMU is operating at a frequency dial is a multiple of that of odier 
subsystems. In this case, it is possible that the IMU is capable of streaming multiple measurements of the 
integrated non-gravitational acceleration and angular velocity over a single system time-step. Consider 
the accumulated measured integrated non-gravitational acceleration as 

:= y>v'„,< 

<= 1 

where it has been assumed diat there are n subdivisions of the interval {tk-i t k ] with a single subinterval 
defined by {U~\ U} and w-here A\^. i £ is given by Eq. (32). Let 

n n n 

Ev£ := >J Avf , bs AV := V b A „ , and ws AVlk := V w A ».i 

<=i fei f=i 
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such that 


£ v m,k = Evjr + b EiV + w EaV|1 . 

Similarly, consider the accumulated measured integrated angular velocity as 

n 

s e c m k ■.= V 

t~ i 

where it has been assumed that there are n subdivisions of the interval [>k-i <*} with a single subinterval 
defined by [t t ] te) and where A0J, is given by Eq. (38). Let 

n n n 

:= byj) := "S' b A e , and wso,k := V w a# / 

<=i *=i <=i 

such tliat 

Ld' n k = l + bse + wst.it- 

3.3 Dynamical Model in Discrete-time 

Define the quaternion A q\(t) representing the rotation from the a priori attitude, that is the attitude 
at time i k -i denoted by 1 > as 

Aq,-(t) := q*(t)®^,k - i 

Let 0 b (t) be the rotation vector parameterization of A q b {t). Assume a small time step in the state prop- 
agation with the 1MU from tk~i to tk- Linearizing Bortz’s equation around small angles and assuming 

“*/«(*)= “‘/(^yields 

0 6 W = < i,k + \o b (t) x u b b/ik tk -1 < t < tk 

which has the solution 

0 k (<)=u>J A ,k(t-fi_i) tk-l <t<t k . (39) 

Letting t -> tk and defining Atk := tk — tk-i yields 

®k = < * , 6/i,k‘^ t <= = ^®k> 

where we assume the IMU possesses accumulated integrating gyros. Therefore, the discretized quater- 
nion propagation is given by 

qt,k = q. 6 (E0t)®9tk- 1 (40) 

where 

co.(ii|Eeti|) J- 

Note that the IMU output T:0l is in the IMU case frame, so when utilizing Eq. (40) we must first transform 
the output using die constant (and known) T“ as T,d h k = T£E0k- As with the attitude discretization, as- 
sume a small time step in the position and velocity propagation with the IMU from tk - 1 to tk. Over this 
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small time step At*. := tk - t*_i, assume that the non-gravitational acceleration a c ng (t) = a.': ng k and the 
gravity a‘ s (t ) = if k _ 1 := g(r i imu j c _ 1 + Tj , A ._ 1 r‘ 9/imu ), where T^_! := Since we assume small 

angles over A tk, we can use Eq. (24) to obtain 

T‘(9?) = (isxs - [o‘x]) < t < t k . (41) 

With T|(<^) = Tj”‘ ($) and = u b b/i k , write Eq. (20) as 

vLAt) = 0L1 + Kk-i Tc«n 5 ,/c - TiU.jITja^xK/^t - tfc— i). (42) 


In Eq. (42), the constant gravity is approximated by g kl x g‘( r| mu t l ). Integrating Eq. (42) twice and 
once, respectively, yields the position and velocity evolution over the interval t*_i < t < tk as 


r imu(9 — T imu,k-1 + v imujc-l(t ~ tk- 1) + ~ **- 1) 2 + ~1'bjc-l^c a ng,k(t ~ f k-lj~ 




v imu (0 - 1 + 9t-l(* _ tk-l) + T b k _ 1 T > ‘ c a^ ig k (t — tk- 1) — ^Tb.k-l [1c o n0,* x ]w»/l,t(* ~ tk- 1) 2 

Letting t it, we have 

r imu,k ~ t'imu^k— 1 T v imu,k-l^tk + ^9fc-lAlfc + ~^’t'),,k l’^c u ng,kAtk — gTjffc— 1 [T7a nff ft x Affc 

= v imu,k-i + 9t-i^ ( tc + Zl* At* - fc-jTjla'g.t xlaij^ t Ai| 

Assuming that the IMU possesses accumulating integrating gyros and accelerometers, we rewrite the 
position and velocity propagation in their final form as 


r 'imu,k ~ r imu,k-l + v \mu,k-l^tk + At^ + ~'^h,k- I 1^3x3 + ~[T c Ed£ X]J T c Et)£Af* (43a) 
vLu,k = vLu,k-i+aL,Atk + Tik -! ( /3X3 + i[T‘E0j x]) T b Evk (43b) 


The discretized mathematical model describing the evolution of the position, velocity, and attitude is 
given by Eq. (40) and Eqs. (43). The estimated states in the EKF are propagated accordingly via 


— t’imu.fe— i +v! mu> k-iAtk + ^g'k-iAtl + - T£ k -i ( hx 3 + -[TjE0£x]') TjEr^Affc (44a) 

= i>Lu,k-i+gLi*tk + Tik-k ( hx3 + i[TjE^x]) T b Ev' k (44b) 

¥i ,k = <?(£#*) ® Qt,k - 1 (44c) 

where g‘ k b := g(r\ mu k _i) and TJ kl := Tj^Lj). Note that the accumulated integrated output of 
the IMU, denoted by E 6 k and Ev k , is corrected with the estimated IMU gyro and accelerometer biases to 
obtain the values Ed k and Et>£ for use in the EKF propagation. 


3.4 State Estimation Errors 

3.4. 1 State Estimation Error Propagation 

Define the multiplicative attitude estimation error at any time t k as 

tq b i,k :=<??,* ®&‘ ■ 
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For small angles, the attitude error is approximately twice the vector part of the quaternion &q h ik , so 

ee.it := ^q\ k . 

Also define 

:= E0£ - £ B b k and e^,,* := bje - 

Using E 0 k to compute 

sin(i||E0‘||)S0‘/||Ee‘|| 


lUm) = 


■ (Mu) 


we denote the associated transformation matrix by T(Y,O h k ). The attitude estimation error propagation 
(from t k - 1 after the previous state update in the EKF to t k prior to the next state update) is given by 

el k = T(E0‘)e+,_ 1 - T*e+ , t _, - Tjwajk. (45) 

Define the position and velocity estimation errors as 

e r ,k '■= fimu , k k and '■= v\ muJe - vj mntk 

Also define 

es v,k '■= Ei>£ - ED* and e^ t k •= bs AV - bs AVj jfc, 

where Eu£ = is transformed from the 1MU case frame. The position and velocity estimation errors 

propagate according to 

e r,k ~ e r,k-l + e t,k-l^k — ^ x e e,fc- 1 (45) 

- + H« w £ AV.» + R 3 v 'X0.k 

e:,k = <*_1 - (tU-i x Efr£)xA e+^ 

- V a et^j,_ 4 + ^ w Ea« V s yfz»,k 

where 


1-, 


1. 


Ha = A-! ( 0*3 + tt[E< X ] T‘Aty , R g = 7f 2j,*_ 1 [E6ix]T*Ai* , 


/ 


1, 


1a, 


v. = ^3X3 + 5 P»txn Tj . and = -T^^lE^x]^. 

The IMU and sensor biases are modeled as random constants. Define the estimation errors in the ve- 
locimeter, attitude camera, and altimeter biases, respectively, as 

e K.„k 


so that 


kyei,* 

and 

e (Ll,k-l : = 

Kc,k 

and 

e K„k- 1 : = 


and 

e b.u,k- 1 

e b Sv ,k 

= 

e tc„Jk- 1 

,k 

= 

e tv.,k- 1 

e b vtl ,k 

= 


e b aet k 

- 

e b a0> k- 1 

e b«u,k 

= 

p + 

e K„Jc- V 
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3.4.2 State Estimation Error Covariance 

Defining the estimation error covariance by P at time tk prior to the state update and P^_ x at time 
tk ~ i after die state update and using the error equations in Eq. (45) and Eq. (46) we obtain the state error 
covariance propagation in Eq. (11) where the state estimation errors are 


e r ,k- 1 

+ 

&r,k 

ev,k-i 


&v,k 

e#,k- 1 


e 0,k 


and el := 


1 






e ^«e.fc-l 





. e *ale,fc . 


and the state transition matrix e R 22x22 and process noise matrix Qk € R 22 * 22 are given by 


hx3 

hx3^k 

-5 (*!U- 1 P»tx]+ l^k-rK^ x Eti‘)x]) A t k 

~Ra 

Rg 

03x3 

03x3 

0 

03x3 

hx3 

~ + jTtU-iP ® l * S«‘)x]) 

-V* 

n 

03x3 

03x3 

0 

03x3 

03x3 

T(S0t) 

03x3 

rpi 

1 C 

03x3 

03x3 

0 

03x3 

03x3 

03x3 

^3x3 

03x3 

03x3 

03x3 

0 

03x3 

03x3 

03x3 

03x3 

hx3 

03x3 

03x3 

0 

03x3 

03x3 

03x3 

03x3 

03x3 

•^3x3 

03x3 

0 

03x3 

03x3 

03x3 

03x3 

03x3 

X 

CO 

O 

h*3 

0 

03x3 

03x3 

03x3 

03x3 

03x3 

03x3 

03x3 

1 


and 

Qk = M k WMl , 

respectively, where 


Mk 


-R a R„ 
-V* V s 
0ix3 -T* 

03x3 03x3 

03x3 03x3 

03x3 03x3 

03x3 03x3 

0 0 


€ R 2 


W = 


W Au 

03x3 


03x3 

w A * 


€ K 6x6 


and 


IQ = ! ( / 3 * 3 + j TjA t k , R g = , 

v« = (l 3x3 + I[E el x]j T b c , and V g = If 6 ‘, x]T* 
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4 Multi-mode Extended Kalman Filter 

4.1 Mathematical model 

The state vector is defined to be 


t'imu.k 

V* i 
tmu.fc 

« 

bs AV 

hnel 

hoc 

bait 


er 


Tiie sensor biases, denoted by b ve i, b QC , and b a i t , for die velocimeter, attitude camera, and altimeter, 
respectively, are modeled as random constants. One caveat to die state estimate update relates to the at- 
titude update. Due to the fact diat we implement a multiplicative update for the attitude quaternion, the 
state vector is R 23 , whereas die state estimation covariance is R 22 * 22 . All of the preceding developments 
have focused on die attitude error in the form of small angles, denoted by 6a. This is done because the 
quaternion of rotation is subject to a unity norm constraint whereas the small angles are not. Therefore, 
if die portion of x£ pertaining to die attitude is given by then the quaternion update is 






(47) 


where q k is the a priori estimate of the quaternion. Eq. (47) yields a quaternion which satisfies the unity 
norm constraint to first order, and so to ensure that the quaternion remains unity norm a brute-force 
normalization is performed such that </,' 1 


4.2 mEKF Hierarchical Structure 

A mEKF hierarchical structure (mEKFHS), shown in Figure 4 is proposed to represent a convenient vehi- 
cle for efficiently collecting data for input and output of the filter. The mEKFHS is comprised of five major 
substructures: meas, imuk, Hoad, Estiinatekm, and Estimatekp. In addition, for simulation purposes, a 
Truth structure is proposed to allow a perfect navigation mode. 

rneask Structure The meask structure passes the measurement data to the mEKF. The structure is com- 
prised of the measurement time tag, f k , the data valid flag, and the measurements. The structure is de- 
tailed in Table 1 . 

imuk Structure The imuk structure passes the IMU Av and AO data to the mEKF. Since the IMU internal 
rates are higher than the mEKF update rate, we allow for the IMU to bin multiple Av and AO values and 
pass them as a group to the mEKF. The suncture is comprised of a variable describing die number of Av 
and AO values being passed and the actual Av and AO values. The structure is detailed in Table 2. 

Truth Structure The Truth structure is the connection to the environment and provides the pathway to 
the true states. In the actual implementation of the mEKF, the Truth structure would not exist. However, 
in the environment framework, it makes sense to have the capability to run a so-called perfect navigation 
mode. The Truth structure provides the necessary link to the environment to allow the mEKF to output a 
state estimate equivalent to the true state. The Trudi structure is detailed in Table 3. 
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Table 1 : The meask structure. 


Sensor 

Variable 

Comment 

Velocimeter 

meask.vr.valid 
meask. vr _ time 
meask. vr 

1 = measurement valid; 0= measurement invalid 
Measurement time tag 

Velocimeter measurement (3-axes in sensor frame) 

Attitude Camera 

meask. ac .valid 
meask. ac .time 
meask. ac 

1 = measurement valid; 0= measurement invalid 

Measurement time tag 

Attitude measurement (quaternion) 

Attitude Camera 

meask. alt. valid 
meask.alt . time 
meask. alt 

1 = measurement valid; 0= measurement invalid 
Measurement time tag 

Altitude measurement (above reference sphere) 

TRN 

meask. tm_ valid 
meask. tm _ time 
meask. tm 

1 = measurement valid; 0= measurement invalid 
Measurement time tag 

TRN measurement of position relative to ground feature 


Table 2: The inuik structure. 


Variable 

Comment 

imumk.nimu 

imumk.deltav 

imumk.deltatheta 

Number of Aw and AS values being passed 
Acs in the 1MU case frame 
ASs in the 1MU case frame 


Table 3: The Truth structure. 


Variable 

Comment 

Trutht 

Simulation time 

Truth.r 

The components of the vehicle position vector in the planet-centered 
inertial coordinate system 

Truthv 

The components of the inertial vehicle velocity vector in the planet-centered 
inertial coordinate system 

Truth, rcg 

The vehicle center of gravity location along the body axes 

Trutha 

Components of measurable (sensed) vehicle acceleration in the planet-centered 
inertial coordinate system 

Truth, g 

The gravity acceleration vector components along the planet-centered inertial 
coordinate system 

Truth.ab 

Tmth IMU accelerations, components of measurable (sensed) acceleration in the 
body coordinate system w/ origin at the IMU 

Truth, w 

Truth IMU angular rates, vehicle angular velocity components in the body 
coordinate system w/ origin at the IMU. 

Truth vr 

The vehicle velocity relative to the rotating planet 

Truth alt 

The vertical altitude of the vehicle above the oblate planet 

Truth thr 

The net thrust (vacuum thrust corrected for atmospheric backpressure effects) 

Truth, la t 

The geocentric latitude of the vehicle. 

Trudr-lng 

The longitude of the vehicle measured east of the prime meridian 

Truthm 

The current vehicle mass 

Imlhq 

Quaternion expressing Ihe rotation from inertial to body 
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Figure 4: mEKF Hierarchical Structure 


Iload Structure The Hoad protides the required mEKF initialization data. Note that the initial state esti- 
mate and initial state error covariance matrix are initialized with the Estimatekm and Estimatekp struc- 
tures. The Iload structure is detailed in Table 4. 
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Table 4: The Iload structure. 


Variable 

Comment 

Iload. nx 

Number of states in the filter 

Iload. np 

Number of states in the covariance 

Iload.nin 

Number of measurements 

Iload. navjate 

Update rate of the mEKF in Hz 

Iload.nav.rdt 

Inverse of the navigation rate in sec 

Iload.perfectnav 

1 = Perfect nav is ON; 0 = mEKF is ON 

Iload. gravity model 

Only choice at current time is ’Central’ 

Iload.planet 

Current choices are ’moon’ or 'mars’ 

Iload. qXs 

Orientation of planet-fixed reference frame with respect 
to die planet surface reference frame 

Iload. r.rei.l 

Location of die surface reference frame origin in the 
planet fixed-frame 

Iload.omegaJLf 

Spin axis of moon at epoch 

Iload.qJeJ 

Orientation of the moon inertial reference frame with respect 
to the planet-fixed reference frame at epoch 

Iload.te 

Epoch associated with qje _f in seconds 

Iload.qimu 

I MU process noise matrix 

Iload.qtune 

liming process noise matrix 
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5 Sensor Modeling 

The mEKF navigation algoritlun is a dual-stale extended Kalman filter (EKE). The EKP is a model- 
based algorithm requiring models of the sensors. In general, there are two classes of sensor models re- 
quired for an analysis of the six degree-of-freedom (6 DOF) integrated guidance, navigation, and control 
(GN&C) system. These are (1) high-fidelity models to support 6 DOF simulation of the sensors, and (2) 
sensor models for the navigation algorithm. This document is concerned only with sensor models for the 
navigation algorithm. To that end, a navigation sensor model includes the following: 

1. a mathematical model represented by a nonlinear equation(s) as a function of die states of die sys- 
tem (such as position, velocity, attitude, etc.). 

2. a measurement mapping matrix comprised of the partial derivatives (of the model above) evaluated 
at the most recent state estimate. 

3. an error model comprised of random noise and systematic errors, including representative values 
of the uncertainty in die various error sources. 

The sensors currently under consideration are: inertial measurement unit (IMU), spherical altimeter, 
velocimeter, and star camera. 

5.1 Sensor hierarchical structure 

A sensor hierarchical structure (SHS), shown in Figure 2 is proposed to represent a convenient vehicle for 
representing the multitude of sensors required for EDL. The SHS is comprised of four major substruc- 
tures: Universe, Nucleus, Identity, and F.rrors. 

Universe Structure The Universe structure provides the parameters required by the filter that describe 
the external environment. At the present time, the Universe structure models Mars and the Moon. For 
Mars we provide the polar radius, equatorial radius, 12, gravitational constant, rotation rate, and the in- 
verse of the principal axes. The variables are described in Table 5. As die models increase in complexity, 
in addition to the parameters currently included, the Universe stiucture might also include, for example, 
any topographic maps, location and orientation of planetary bodies (such as die sun), and transforma- 
tions from planet fixed to inertial reference frames. 

Nucleus Structure The Nucleus structure describes die key elements of die sensor around which die 
model is constructed including, for example, the rate (internal), output rate, on/off events, and spe- 
cial processing (hardware specific) parameters. Currentiy, the Nucleus structure models die velocimeter, 
spherical altimeter, attitude camera, and IMU. Nominal values for each parameter is si '.nun in Table 6. 

Identity Structure The Identity structure describes die physical aspects of the sensor such as, for exam- 
ple, the spacecraft sensor layout (sensor locations, sensor reference frame orientations, etc.), and beam 
locations (for the velocimeter and altimeter). 

F.rrors Structure The Errors structure includes the random and systematic errors associated with each 
sensor, such as biases, noise slrengdis, misalignments, and time constants. 
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lable 5: Universe structure. 


Name 

Comment 

Universe.mars.Ae 
Universe. marsApo 
Universe, mars .Aeq 
Universe. mars.J2 
Universe, mars . Mu 
Universe, mars. Omega 
Universe, mars.prcpllengths 

Radius of Mars for gravity calculations 
Polar radius of Mars 
Equatorial radius of Mars 
Gravity J2 term of Mars 
Gravitational constant of Mars 
Rotation rate of Mars 
Computed value 

Universe, moon. Ae 
Universe, moon. Apo 
Universe. moon.Aeq 
Universe, moon. J2 
Universe, moon. Mu 
Universe, moon. Omega 
Universe, moon. prcpUengths 

Radius of the Moon for gravity calculations 
Polar radius of the Moon 
Equatorial radius of the Moon 
Gravity J2 term of the Moon 
Gravitational constant of the Moon 
Rotation rate of the Moon 
Computed value 


Table 6: Nucleus structure. 


Name 

Comment 

Sensor.Nucleus.velocimeter.edit-onoff 
Sensor. Nudeus.velocimeter.edit_gate 
Sensor.Nucleus.velocimeter.underweight_onoff 
Sensor.Nudeus.velocimeter.underweight_accept 
Sensor.Nucleus.velocimeter.rate 

Velocimeter residual edit enabled/disabled 
Specification of probability gate (e.g., 9.14 * 99.75%) 
Velocimeter underweighting enabled/disabled (1/0) 
Underweighting scale factor 0 < p < 1 
Velocimeter measurement rate in Hz 

Sensor.Nucleus.altimeter.spherical.edit-onoff 
Sensor. Nucleus.altimeter.spherical.edit.gate 
Sensor. Nucleus. altimete r. sphe rical.undenveight_onoff 
Sensor. Nucleus. altimeter, spile ri cal. unrienveight_acrep[ 
Sensor.Nucleus.altimeter.spherical.rate 

Altimeter residual edit enabled/disabled (1/0) 
Specification of probability gate 
Altimeter underweighting enabled/disabled (1/0) 
Underweighting scale factor 0 < p < 1 
Alfimeter measurement rate in Hz 

Sensor. Nucleus. attitudecamera.edtt_onof( 
Sensor.Nucleus.attitudecamera.edit_gate 
Sensor.Nucleus.attitudecamera.underweight-onoff 
Sensor. Nucleus. attitudecarnera. u tide rweighLaccept 
Sensor.Nucleus.attimdecamera.rate 

Star camera residual edit enabled/disabled (1/0) 

Specification of probability gate 

Star camera underweighting enabled /disabled (1/0) 

Underweighting scale factor 0 < p < 1 

Attitude measurement rate in Hz 

Sensor.Nucleus.imu.deltav.threshJimit 

Sensor.Nucleus.imu.deltav.thresh-onoff 

Sensor.Nucleus.imu.deltav.thresh_sigsq 

Sensor.Nucleus.imu.deltatheta.threshJimit 

Sensor.Nucleus.imu.deltatheta.thresh_onoff 

Sensor.Nucleus.imu.deltatheta.thresh_sigsq 

Sensor. Nucleus.imu.rate 

Sensor. Nucleus.imu.nimu 

Tolerance for AV r thresholding 

IMU AV thresholding enabled /disabled (1/0) 

Standard deviation of expected IMU AE noise 

Tolerance for A0 thresholding 

IMU A# thresholding enabled /disabled (1/0) 

Standard deviation of expected IMU A0 noise 

IMU internal rate in Hz 

Number of buffered IMU data 
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Table 7: Identity Structure. 


Name 

Comment 

Sensor.Identity.altimeter.spherical.r.altimu.b 

Sensor.Identity.velocimeter.r_velimu_b 

Sensor. Identity.:! itiuidecameia.rjrimuJ) 

Sensor.Identity.q_c_b 

Sensor. Identity.q-ti-b 

Sensor. Identity.qjtcJb 

Sensor.Identity.q.vr.b 

Location of the altimeter wrt die I MU 
represented in the body reference frame 
Location of the velocimeter wrt the IMU 
represented in the body reference frame 
Location of the attitude camera wrt the I MU 
represented in the body reference frame 
Attitude quaternion of IMU case frame to the 
body reference frame 
Attitude quaternion of altimeter reference 
frame to the body reference frame 
Attitude quaternion of attitude camera 
reference frame to the body reference frame 
Attitude quaternion of velocimeter 
reference frame to the body reference frame 


Table 8: Errors Structure. 


Name 

Comment 

Sensor.Errors.velocimeter.bias 
Sensor.Errors. velocimeter. rvel 

Velocimeter measurement bias 
Velocimeter measurement noise 
covariance matrix 

Sensor.Eirors.altimeter.spherical.bias 
Sensor.Errors. altimeter, spherical. ralt 

Spherical altimeter measurement bias 
Spherical altimeter measurement noise 
covariance 

Sensor.Errors.attitudecamera.bias 
Sen sor. Error s. a t ti tu de ca me ra . rac 

Attitude camera measurement bias 
Attitude camera measurement noise 
covariance matrix 

Sensor.Errors. imu.accel_spec.bias 
Sensor.Errors. imu.acceLspec_noise 
Sensor.Errors.imu.gyro.spec.bias 
Sensor.Errors. imu.gyro_spec_noise 

IMU Av bias standard deviation 
IMU Av noise covariance 
IMU A 9 bias standard deviation (1 -<t) 
IMU A8 noise covariance 
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For all four sensors, the data is nominally provided every At seconds computed via 


Sensor. Nucleus.velocimeter.specjdt 
Sensor.Nucleus.altimeter.spherical.specjdt 
Sensor . Nucleus, at ti t udecaraera. spec_d t 
Sensor.Nucleus.imu.spec_dt 


1.0/Sensor.Nucleus.velocimeter.rate 
1 .0/Sensor. N ucleus.altimeter. spherical, rate 
1. 0/Sensor. Nucleus.attitudecamera.rate 
1 .0/Sensor. N ucleus.imu.rate 


Also, we have 

Universe.mars. prcpllengths= 
Universe.moon.prcpllengths- 


U niverse.mars. Aeq 
0 
0 

Universe, moon. Aeq 
0 
0 


Universe, mars. Aeq 0 

0 Universe.mars. Apo 

0 0 

Uni verse, moon. Aeq 0 

0 Universe, moon. Apo 


5.2 Inertial Measurement Unit 

With the sensor data structures defined as in the previous sections, we can represent the I MU with the 
Universe, Nucleus, Identity, and Error data structures. Recall that the model of the I MU was presented in 
Section 3.2. 

5.2. 1 IMU Sensor Hierarchical Structure 

The IMU sensor hierarchical structure is shown in Figure 5. In Figure 5, we show the planet described 
in the Universe structure as Mars. Other planets can be represented. For example, the variable Uni- 
verse.mars.Ae would be Universe.moon.Ae for the Moon. 



Figure 5: IMU Sensor Hierarchical Structure 
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5.3 Altimeter 

5.3.1 Mathematical model 

The estimated spherical altimeter measurement at time tk is 

hk = — r sp/ij + b altk , 

where 

"T T b j e r alt / imu 


r nit. k — 


and 


ti;* = t t (^;) 


The value of r sp h is currently input as a constant in the spherical altimeter model, but it could be varied 
to capture the local topography, if desired and if a topography model were available to navigation. If 
an altitude measurement is given at time tk by hk, and the estimate is given by hk, dten the spherical 
altimeter residual is 

r* = hk — hk- 

The measurement sensitivity for the spherical altimeter has the form 

H* = [ III 0ix3 H 3 0jx3 0 [x3 0ix 3 0ix3 Hs ] . 

The specific elements of the measurement sensitivity are 

1 .i~.T 


III 


H 3 = -T 


‘oil,* 

1 


\ l alt,k\ 

H« = 1 

where the estimated position of tire altimeter in the inertial frame is 

fait, A: ^imu,k "b I'fc.&l’atr/inH** 

5.3.2 Altimeter Sensor Hierarchical Structure 

The altimeter sensor hierarchical structure is shown in Figure 6. 


5.4 Velocimeter 

5.4.1 Mathematical model 

The estimated relative velocity measurement is 

^rel = TJTf,* X i' V el,k] + Kel,k’ 

where u>‘ L is the rotation vector of the planet and the estimated position and velocity of the velocimeter 
in the inertial frame are, respectively, 

^vel,k ^imu,k ”b '^- byk^vel/imu 

^vel,k v imu,fc “b * ^vel/imu\‘ 
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Figure 6: Altimeter Sensor Hierarchical Structure 


With the relative velocity measurement given at time tk by v” el k , and the estimate given by v” el k , then 
the velocimeter residual is 

T k = V rel,k ~ Vreljc- 

The measurement sensitivity for the relative velocimeter has the form 

H/c = [ Hi Ho H3 03x3 03 x 3 He O3X3 03x1 ] . 

The specific elements of the measurement sensitivity are 

Hi = -TUTQu^x] 

H, = TJTf* 

H3 = -TX xr‘ ei/im Jx] +^tr t [ w ix]Tt; ic [r^ /i m U x] 

+ T t ['it.fc - [ u 'l x fvel.fel] x ] 

H 6 = I 3 X 3 

5.4.2 Velocimeter Sensor Hierarchical Structure 

Tite velocimeter sensor hierarchical structure is shown in Figure 7. 

5.5 Attiude Camera 

5.5.1 Mathematical model 

The estimated quaternion attitude camera measurement is 

kVv = q hjuk 0 <JT ® k'tk ® 
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Figure 7: Velocimeter Sensor Hierarchical Structure 


where 


/ ^ ac,k ^qc,A 

¥) 


Qb,Tj,k 


and b~ c k = ||b~ c k \\, where b~ ck is the estimate of die attitude camera bias described in Section'3.4.1. 
Therefore, if a measurement of die quaternion attitude camera is given at time t-k by q°r,/t> and die estimate 
is given by ij“^ k , then the attitude camera residual is 

ffc “ 9sr,tfFlsr,fc "b (Isr.k * f t,,- fc — ttsr.iQsr.fc) ■ 

Ilie measurement sensitivity for a quaternion attitude camera measurement has the form 

H lc — [ 0 3x 3 03x3 H s 03 x3 0 3 x3 0 3 x3 H 7 0 3 xl ] . 

The specific elements of the measurement sensitivity are 


H 3 

h 7 


l b,n,k l b 

I 3 x3, 


where T k ^ is the transformation matrix created by the estimated bias-noise quaternion, 

Tb,,,* = T (Qb,r,,k). 

5.5.2 Attitiide Camera Sensor Hierarchical Structure 

The attitude camera sensor hierarchical structure is shown in Figure 8. 
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Figure 8: Attitude Camera Sensor Hierarchical Structure 
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